Project 5 This assignment can be worked on and submitted in groups of up to three persons. Goals. In this project, we will work with parabolic and hyperbolic partial differential equations, and...

1 answer below »
Number of pages will vary based on length of full report including labeled pictures/diagrams.


Project 5 This assignment can be worked on and submitted in groups of up to three persons. Goals. In this project, we will work with parabolic and hyperbolic partial differential equations, and combine elementary techniques from 2p-BVPs with time-stepping methods from IVPs. The goal is to gain experience with the method of lines and stability constraints of Courant-Friedrichs-Lewy (CFL) type on the time-step ∆t. In the first part we shall primarily work with the parabolic diffusion equation ut = uxx. In the second part we will work with conservation laws such as the hyperbolic advection equation ut = −ux. The objective is to gain an elementary understanding of and experience with “evolution equations” (time dependent partial differential equations), and the many varying properties one encounters in such problems. In addition, we will work with different types of boundary conditions, and it is important to implement boundary conditions correctly in order to get the correct results. We will only work in one space dimension, and on equidistant grids. We will use explicit and implicit time-stepping methods, in order to study the properties of some important elementary finite difference methods. We will also observe some undesirable phenomena such as numerical instability and numerical damping, and learn how to cope with them. Part 1. The diffusion equation The diffusion equation is ut = uxx u(t, 0) = u(t, 1) = 0 u(0, x) = g(x) and its method-of-lines semi-discretization is u̇ = T∆xu 1 on an equidistant grid on [0, 1], with ∆x = 1/(N + 1) and initial condition g(x) sampled on the grid. Generate the grid using linspace, so that you have N internal points on (0, 1). The matrix T∆x is the usual Toeplitz symmetric tridiagonal matrix that approximates ∂2/∂x2 to second order accuracy on the grid. In Project 3 you have generated this (and similar operators. Go back to your old programs and retrieve the code segments you need in order to create T∆x. You will need (some representation of) the operator T∆x in order to solve equation systems involving this matrix. Next, write your Matlab file eulerstep and make sure that it takes one time step, and that it works in the matrix–vector case. It will probably be necessary to make some minor modifications so that you can call eulerstep with a command of the type unew = eulerstep(Tdx,uold,dt) in order to solve u̇ = T∆xu by the Euler method, um+1 = um + ∆t · T∆xum. Task 1.1 (1p) Write a script that constructs T∆x, assuming N interior points on the grid and solves the semidiscretization u̇ = T∆xu by the Euler method. 1. Visualize the computational procedure by plotting the numerical solu- tion u(m∆t, n∆x) over the t, x plane. Basically, you should apply the following procedure. You determine your spatial grid size ∆x and time step ∆t from the number N of internal mesh points on the x-interval [0, 1] and the number M of time steps to go from t = 0 to t = tend, respectively. Always run your experiments by varying N and M , not by varying ∆x and ∆t directly. Then construct a vector xx = linspace(0,1,N+2) and a vector tt = linspace(0,tend,M+1) and run [T,X]=meshgrid(tt,xx) to generate the independent vari- ables for 3D graphics in Matlab. If you are not familiar with these commands, check Matlab’s help. When you solve the PDE, save the solution vectors over the meshed grid. Don’t forget to insert the boundary values! When the computa- tion is done, make a 3D color graph using either the mesh or the surf command. (Run help to find out how they work.) 2 2. Determine experimentally the CFL condition on the time step ∆t for the explicit Euler method applied to the semi-discretization. Choose your initial condition g(x) as you like, but make sure that they are com- patible with the boundary conditions at x = 0 and x = 1, respectively, and make sure that g(x) is not an eigenfunction of ∂2/∂x2. 3. Try to obtain representative plots of the solution both when you fulfill the CFL condition and when you violate it. In the latter case, don’t overdo it; try to see what happens when you violate the CFL condition by a small amount only, so that the numerical instability is clearly visible. Task 1.2 (1p) Next implement the Crank–Nicolson method. This is the same as applying the trapezoidal rule to the equation, un+1 = un + ∆t 2 (T∆xu n + T∆xu n+1). Implement a function unew = TRstep(Tdx,uold,dt) Solve the same diffusion equation as in Task 1. Verify that you can now obtain stable solutions where you previously observed instability due to a violated CFL condition. Try to verify that you can use much larger time steps ∆t than before, but remember to vary time step and Courant number by varying N and M . Check the the results are the same no matter what initial conditions you use. Part 2. The advection equation The linear advection equation is ut + aux = 0, with appropriate initial and boundary conditions. We will only work with periodic boundary conditions. This means that the solution satisfies u(t, 0) = u(t, 1) for all times t. That is, the solution “wraps around” from x = 1 to x = 0, and vice versa. With periodic boundary conditions, we represent the solution u on the entire grid, including the boundaries throughout the entire computation. Let the solution vector at time t be an N -vector and introduce an equidistant 3 grid on [0, 1], such that 0 = x1, . . . , xj = (j−1)∆x, . . . , xN = (N−1)∆x and xN+1 = 1. This implies that ∆x = 1/N , and that we are going to compute the solution at the N points x1, . . . xN . The solution at the right boundary, xN+1 = 1, is uN+1 and is identified with u1 at the left boundary, x1 = 0. Make sure you have understood exactly what this implies (make a sketch of the grid!), as your code will differ from what you did in the parabolic case. Initial conditions are always needed, and have the form u(0, x) = g(x). Note that the initial condition must satisfy g(0) = g(1) and g′(0) = g′(1) when periodic boundary conditions are used. Do not use a trigonometric initial condition such as sin 3πx. One of the simplest methods of order 2 for the advection equation is the Lax–Wendroff scheme. It can be derived using a Taylor series expansion. Thus, u(t+ ∆t, x) = u(t, x) + ∆tut + ∆t2 2! utt + O(∆t 3). From the differential equation ut + aux we get utt = −auxt, as well as utx = −auxx. Therefore we have utt = a2uxx provided that the solution u is sufficiently differentiable. Inserting this into the Taylor series expansion we get u(t+ ∆t, x) = u(t, x)− a∆tux + a2∆t2 2! uxx + O(∆t 3). This gives the Lax–Wendroff scheme, un+1j = aµ 2 (1 + aµ)unj−1 + (1− a2µ2)unj − aµ 2 (1− aµ)unj+1. Note that the method coefficients are not symmetric – instead, they change with the flow direction as determined by the sign of a. Task 2.1 (1p) Implement a Lax–Wendroff solver for the scalar advection equation with periodic boundary conditions. Write a function function unew = LaxWen(u, amu) that takes a step of size ∆t starting from the “initial condition” u. The parameter amu is the product a∆t/∆x. Write a script for running your Lax– Wendroff solver. Test your solver using some pulse-like initial data, e.g. g(x) = e−100(x−0.5) 2 . 4 Check how this pulse propagates for 5 units of time, and verify that the method works both for positive and negative a. Pay close attention to how you implement the periodic boundary conditions (sketch the grid and the computational stencil before you write the program). Produce plots of the solution for at least two different CFL numbers. In addition, plot the L2 norm (RMS norm) of the solution vs time t ∈ [0, 5] for a∆t/∆x = 1 as well as for a∆t/∆x = 0.9, and explain the difference. Is the Lax-Wendroff method conserving the norm of the solution? Does the amplitude remain constant, or is it damped? Motivate why one actually wants to run the code at the CFL stability limit. 5 Matlab Script Codes Prj5p1 clear all;close all; %Part 1a: The Diffusion Equation % % ut= uxx % u(t,0)= u(t,1)= 0 % u(0,x)= g(x) % % dx= 1/(N+1) % u_t= Tdx*u ~ u^(m+1)= u^m + delta_t*Tdx*u^m <--- this="" is="" the="" explicit="" euler="" form="" %="" %="" exact="" solution="" ga="@(x)" sin(pi*x)="" -="" sin(3*pi*x);="" ua="@(x,t)" exp(-(pi)^2*t)*sin(pi*x)="" -="" exp(-(3*pi)^2*t)*sin(3*pi*x);="" %="" spatial="" finite="" difference="" discretization="" n="30;" dx="1/(N+1);" %="" equidistance="" grid="" on="" [0,1]="" with="" n+1="" subintervals="" t="diag(-2*ones(N,1),0)" +="" diag(ones(n-1,1),-1="" )+="" diag(ones(n-1,1),1);="" tdx="1/(dx^2)*T;" xx="linspace(0,1,N+2);" %="" equidistance="" grid="" on="" [0,1]="" with="" n+1="" points="" xx="xx';" %="" temporal="" euler="" explicit="" tend="0.1;" mg="[200," 50];="" for="" i="1:length(Mg)" m="Mg(i)" dt="tend/M" %="" total="" m="" time="" stepping="" tt="linspace(0,tend,M+1);" %="" time="" instants="" z="zeros(N+2,M+1);" %="" soution="" u(x_i,t_m)="" za="zeros(N+2,M+1);" %="" exact="" solution="" err="zeros(M+1,1);" z(:,1)="ga(xx);" %="" initial="" condition="" u(x_i,0)="g(x_i)" za(:,1)="ua(xx,tt(1));" err(1)="norm(Za(:,1)-Z(:,1));" for="" m="1:M" z(:,m+1)="eulerstep(Tdx,Z(:,m),dt);" %="" explicit="" euler="" za(:,m+1)="ua(xx,tt(m+1));" err(m+1)="norm(Za(:,m+1)-Z(:,m+1));" end="" figure(i)="" [t,x]="meshgrid(tt,xx);" subplot(3,1,1)="" surf(t,x,za)="" %="" this="" is="" the="" 3d="" plot="" of="" the="" exact="" solution="" subplot(3,1,2)="" plot(tt,err,'o-')="" %="" this="" is="" the="" 2d="" plot="" of="" the="" errors="" between="" the="" exact="" solution="" and="" computed="" ones.="" subplot(3,1,3)="" surf(t,x,z)="" %this="" is="" going="" to="" give="" us="" the="" 3d="" plot="" cfl="dt/dx^2" nrm0="norm(Z(:,M));" nrm1="norm(Z(:,M+1));" if="" (cfl="">< 1/2)="" title('cfl="" condition="" are="" met="" and="" solutions="" should="" be="" stable')="" elseif="" (cfl="">1/2 & nrm11/2 & nrm1>nrm0) title('CFL Condition are not met and solutions are not stable') end end Prj5p1b clear all;close all; %Part 1.2: Solving the Diffusion Equation by the Crank-Nicolson Method % % ut= uxx % u(t,0)= u(t,1)= 0 % u(0,x)= g(x) % % Spatial Finite Difference Discretization
Answered 2 days AfterApr 29, 2022

Answer To: Project 5 This assignment can be worked on and submitted in groups of up to three persons. Goals. In...

Pawan answered on May 02 2022
101 Votes
Objective
In this project, we will work with parabolic and hyperbolic partial differential equations, and combine elementary techniques from 2p-BVPs with time-stepping methods from IVPs. The goal is to gain experience with the method of lines and stability constraints of Courant-Friedrichs-L
ewy (CFL) type on the time-step ∆t.
The diffusion equation
In this part, we solve the 1D diffusion equation given by
with the boundary conditions given by
.
For the current project, we choose
We solve the above partial differential equation (PDE) using the method-of-lines. In this method the given PDE is descretized as
,
where is the Topelitz matrix which approximates the second derivative of the function . We integrate the above equation using two different integrators and examine its behaviour.
For both the cases, we consider equidistant grid points.
Euler Step
We integrate the discrete equation using
We use the grid point N = 30 and simulate the problem for two time steps. We evaluate the time step by dividing the final time i.e. 0.1 sec with 50 and 200. In Figure 1, we plot the result at t=0.1 sec for .

Figure 1: The exact solution (top), error with time (middle) and the computed solution (bottom) for the 1D diffusion problem for .
The top and bottom plot shows the value of the wave at t=0 to t=0.1 sec for the entire domain of size [0, 1]. The z axis shows the magnitude of the wave at any time t and location x.
Clearly, the solution is unstable for . In order to make the solution stable, we use a small timestep of size . In the Figure 2, we plot the exact solution, the evolution of error and the computed solution at t=0.1 sec.
Figure 2: The exact solution (top), error with time (middle) and the computed solution (bottom) for the 1D diffusion problem for .
Clearly the solution is stable till 0.1 sec. The error remains bounded. Therefore, it is recommended to used small timestep in order to satisfy the CFL criteria.
Crank–Nicolson method
In this section, we change the integrator applied in the last problem, We use an explicit integrator. The integration scheme is given by
    
One of the benefits of using an explicit integrator is that it is not limited by stability criteria. In order to show the benefit of the above scheme, we use a time step of 0.1/40=0.025 sec. This time step is even higher than the previous highest case.
Figure 3: The exact solution (top), error with time (middle) and the computed solution (bottom) for the 1D diffusion problem for using the Crank-Nicolson integrator.
In figure 3, we plot the exact solution, the error evolution with time and the numerically commuted solution. Clearly the solution is accurate even for a very coarse time step. Therefore, we can conclude that explicit integrators are not limited by any stability criterion.
                           
            
                
                    
The advection equation
In this section, we solve the advection equation given...
SOLUTION.PDF

Answer To This Question Is Available To Download

Related Questions & Answers

More Questions »

Submit New Assignment

Copy and Paste Your Assignment Here