Lax-Wendroff and Crank-Nicolson Matlab simulation of linear advection

Question

I am implementing two schemes: Lax-Wendroff and Crank Nicolson to solve the PDE linear advection in Matlab. My simulation seems to run; however, between the transition in time space the numerical solutions obtained do not seem to look right.

This is the schemes that I implement:

In Matlab implementation, I am trying to solve $v_{j}^{n+1}$ which is denoted as v_new. For Lax-Wendroff, I implement it as the for-loop updating and I take into account the periodic solution. For Crank-Nicolson, I make a matrix for both $v_{j}^{n+1}$ and $v_{j}^{n+1}$ because it is implicit method.

The following is my Matlab implementation. I cannot fix as I describe above.

% the advection speed
a=1.;
% number of grid points in space (not counting the point at x=1, which is the same as the point at x=0, due to the periodicity)
N=51;
% spatial step
h =(1./(N));
% safety constant for stability (should be smaller than 1)
cfl=.8;
% CFl timestep limit for the explicit methods
dt=cfl*h/a;
% time at which we want to end the simulation
t_end=1.;    
% number of timesteps to be taken
n_it=t_end/dt;
% number of methods: Lax-Wendroff vs Crank-Nicolson. 
n_methods=2;
% initialize some arrays
% spatial grid (not including the point at x=1, which is the same as the point at x=0, due to the periodicity)
x=(0:h:1-h);
% temporal grid
t=(0:dt:n_it*dt);
% arrays for the numerical approximations at times n and n+1
v_new=zeros(N,n_methods);
v_old=zeros(N,n_methods);
% array for the exact solution
v_exact=zeros(N,1);
    
% the initial condition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n_methods
    v_old(:,i)=exp(-50*(x-.5).^2);
end
[ACN,BCN]=mat_linadv_CN(N,a,h,dt);  




    
% the main iteration loop
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
for iter = 1:n_it
% method 1: Lax-Wendroff(treat the first and last separately into account the periodicity)
for i=2:N-1
    v_new(i,1)=v_old(i,1) + (((a^2)*dt)/2*h^2)*(v_old(i+1,1)-2*v_old(i,1)+ ... 
                         v_old(i-1,1))-(a*dt)*(v_old(i+1,1)-v_old(i-1,1));
end 
    v_new(1,1)=v_old(1,1) + (((a^2)*dt)/2*h^2)*(v_old(2,1)-2*v_old(1,1)+...
                    v_old(N,1))-(a*dt)*(v_old(2,1)-v_old(N,1)); 
        v_new(N,1)=v_old(N,1) + (((a^2)*dt)/2*h^2)*(v_old(1,1)-2*v_old(N,1)+...
                   v_old(N-1,1))-(a*dt)*(v_old(1,1)-v_old(N-1,1));
% method 2: Crank-Nicolson. 
v_new(:,2)=ACN\(BCN*v_old(:,2));
% the exact solution with a=1
v_exact(:) = exp(-50*(x-a*t(iter+1)-.5).^2) + ...
             exp(-50*(x-a*t(iter+1)+.5).^2);
% graphical output
plot(x,v_exact(:),'^r-')
hold on

for i=1:n_methods
    plot(x,v_new(:,1),'*b-')
    plot(x,v_new(:,2),'+g-')
end
axis([0 1 -.2 1.2])
xlabel('x')
ylabel('v')
title('linadv simulation')
legend('exact','Lax-Wendroff','Crank-Nicolson'); 
grid on;
grid minor; 
hold off
pause(0.001)

% prepare for the next iteration

v_old=v_new;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%implicit method Crank-Nicolson (CN): 
function [ACN,BCN]=mat_linadv_CN(N,a,h,dt)
e=ones(N,1); 

const1=(-a/(4*h));
const2=(1/dt);
const3=(a/(4*h));

ACN=spdiags([const1*e const2*e const3*e],[-1 0 1], N,N); 
ACN(1,N)=const3;
ACN(N,1)=const1; 

BCN=spdiags([const1*e const2*e const3*e],[-1 0 1],N,N); 
BCN(1,N)=const3;
BCN(N,1)=const1;
end

Answer 1

I have fixed the Lax-Wendroff part, see code below. There were several typos in the coefficients and in the expression of the scheme, see Wikipedia. These remarks might help you fix the whole thing. By the way, note that alternative Matlab implementations are available on this site, see related posts.

% the advection speed
a=1.;
% number of grid points in space (not counting the point at x=1, which is the same as the point at x=0, due to the periodicity)
N=51;
% spatial step
h =(1./(N));
% safety constant for stability (should be smaller than 1)
cfl=.8;
% CFl timestep limit for the explicit methods
dt=cfl*h/a;
% time at which we want to end the simulation
t_end=1.;    
% number of timesteps to be taken
n_it=floor(t_end/dt);
% initialize some arrays
% spatial grid (not including the point at x=1, which is the same as the point at x=0, due to the periodicity)
x=(0:h:1-h);
% temporal grid
t=(0:dt:n_it*dt);
% arrays for the numerical approximations at times n and n+1
v_new=zeros(N,1);
v_old=zeros(N,1);
% array for the exact solution
v_exact=zeros(N,1);

% the initial condition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v_old(:)=exp(-50*(x-.5).^2);

% the main iteration loop
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
for iter = 1:n_it
% method 1: Lax-Wendroff(treat the first and last separately into account the periodicity)
for i=2:N-1
    v_new(i)=v_old(i) + (a^2*dt^2/(2*h^2))*(v_old(i+1)-2*v_old(i)+ ... 
                     v_old(i-1))-(a*dt/(2*h))*(v_old(i+1)-v_old(i-1));
end 
    v_new(1)=v_old(1) + (a^2*dt^2/(2*h^2))*(v_old(2)-2*v_old(1)+...
                    v_old(N))-(a*dt/(2*h))*(v_old(2,1)-v_old(N)); 
    v_new(N)=v_old(N) + (a^2*dt^2/(2*h^2))*(v_old(1)-2*v_old(N)+...
                    v_old(N-1))-(a*dt/(2*h))*(v_old(1)-v_old(N-1));
% the exact solution with a=1
v_exact(:) = exp(-50*(x-a*t(iter+1)-.5).^2) + ...
             exp(-50*(x-a*t(iter+1)+.5).^2);
% graphical output
plot(x,v_exact(:),'^r-')
hold on

plot(x,v_new(:),'*b-')
axis([0 1 -.2 1.2])
xlabel('x')
ylabel('v')
title('linadv simulation')
legend('exact','Lax-Wendroff'); 
grid on;
grid minor; 
hold off
pause(0.001)

% prepare for the next iteration

v_old=v_new;

end