Having Issues With Discretizing Lax Wendroff Scheme

Hello,

I am trying to model a concentration of some substance using 1D Advetion and Diffusion. I am comparing ther Upwinding Scheme and the Lax Wendroff Scheme to the analytic solution for the concentration. Both schemes use the same initial conditions and have the same domains. My issue is when I try to discretize the Lax-Wendroff Scheme I don't get an array back but a single value. This is a problem for my code since it uses periodic boundaries and it breaks when there isn't an array input in the bounds. I am not exactly sure where the problem comes from but I believe it is from my discretization equation for the Lax Wendroff scheme.

The original equation I used for my discretization is :

where C is the concentration, x is the position in the x direction, t is time, D is the diffusion coefficent, u is the advection velocity, the i subscripts keep track of the position in the x direction, the n superscripts keep track of the position in time.

I rearranged it to solve for C(i)^(n+1):

The inital condtion for this scheme is:

and the analytic solution is:

For this code the input values are D = 0.001, u = 0.5, and L =2

Here is my code and the problem is occuring at my Periodic Bounds of Clwnew(1) = Clwnew(Nx-1), this is because the Clwnew value is coming up as a single value instead of an array and I am not sure what I need to do to get to come up as an array after its discretization. Any help would be greatly appreicated.

%% 1D Advection - Diffusion Problem

% Inputs

L=2; % length of domain

D = 0.001; % diffusion coefficient

Nx=50; % points along x

Nt=1000; % length of time

dt=0.01; % time step

% Create grid

x=linspace(0,L,Nx);

dx=x(2)-x(1);

% Initial condition

t=0;

Ci = @(x) cos(2*pi*x/L); % Inital Concentration

CU = Ci(x); % Upwinding

Clw = Ci(x); % Lax-Wendroff

% Define Advection constant

Ux = 0.5;

CUnew = CU;

Clwnew = Clw;

for n=1:Nt

% Update time

t=t+dt;

% Update concentration on interior points

for i=2:Nx-1

% Upwind dT/dx

if Ux>0

dCUdx=(CU(i)-CU(i-1))/dx; % Backward

d2CUdx2 = (CU(i-1)-2*CU(i)+CU(i+1))/(dx^2);

else

dCUdx= (CU(i+1)-CU(i))/dx; % Forwards

d2CUdx2 = (CU(i-1)-2*CU(i)+CU(i+1))/(dx^2);

end

% Update concentration

CUnew(i)=CU(i)+dt*(-Ux*dCUdx + D*d2CUdx2); % UpWinding

Clwnew = Clw(i) + dt*(-Ux*Clw(i+1)-Clw(i-1)/(2*dx) +...

(Ux^2*dt/2)*((Clw(i+1)-2*Clw(i)+Clw(i-1))/(dx^2)) +...

D*(Clw(i-1)-Clw(i)+Clw(i+1))/(dx^2)); % Lax-Wendroff

C_analytic = @(x,t) cos(2*pi*(x(i)-Ux*t(n))/L)*exp(-D(2*pi/L)^2*t(n)); % Analytic Solution

end

% Periodic Bounds

CUnew(1) = CUnew(Nx-1);

CUnew(Nx) = CUnew(2);

Clwnew(1) = Clwnew(Nx-1);

Clwnew(Nx) = Clwnew(2);

% Tranfer new concentration into C array

CU=CUnew;

CLw=Clwnew;

if rem(n,100)>0

figure(2);clf(2);

plot(x,CU)

hold on

plot(x,Clw)

plot(x,C_analytic)

xlabel('x')

ylabel('T(x)')

title(['Time = ',num2str(t)])

legend('Upwinding','Lax Wendroff','Analytic')

drawnow

end

end

[Torsten]

To impose periodic boundary conditions, the loop must run over all points from 1 to Nx and you must set

Clwnew_null = u(nx-1);
Clwnew_nxp1 = u(2);

for artificial points outside the interval of integration.

The way the problem is programmed above only solves in dx:L-dx, not in 0:L.

And

Clwnew = Clw(i) + dt*(-Ux*Clw(i+1)-Clw(i-1)/(2*dx) +...
    (Ux^2*dt/2)*((Clw(i+1)-2*Clw(i)+Clw(i-1))/(dx^2)) +...
    D*(Clw(i-1)-Clw(i)+Clw(i+1))/(dx^2)); % Lax-Wendroff

should read

Clwnew(i) = Clw(i) + dt*(-Ux*(Clw(i+1)-Clw(i-1))/(2*dx) +...
    Ux^2*dt/2*(Clw(i+1)-2*Clw(i)+Clw(i-1))/dx^2 +...
    D*(Clw(i-1)-Clw(i)+Clw(i+1))/dx^2); % Lax-Wendroff

And why do you separate the two terms

Ux^2*dt/2*(Clw(i+1)-2*Clw(i)+Clw(i-1))/dx^2

and

D*(Clw(i-1)-Clw(i)+Clw(i+1))/dx^2

?

Why not

Clwnew(i) = Clw(i) + dt*(-Ux*(Clw(i+1)-Clw(i-1))/(2*dx) +...
    (Ux^2*dt/2+D)*(Clw(i+1)-2*Clw(i)+Clw(i-1))/dx^2)

[Logan Parrish]

Thanks for the response, I do have one question about the periodic boundaries you wrote down, whar is the u value that you are indexing? Also to answer your question about why not combining the terms, I felt it would be a bit less complicated to input as two seperate terms instead of one large one.

[Torsten]

Clw_null = Clw(nx-1);
i=1;
Clwnew(i) = Clw(i) + dt*(-Ux*(Clw(i+1)-Clw_null)/(2*dx) +...
            Ux^2*dt/2*(Clw(i+1)-2*Clw(i)+Clw_null)/dx^2 +...
            D*(Clw(i+1)-2*Clw(i)+Clw_null)/dx^2); % Lax-Wendroff

 for i=2:Nx-1
   Clwnew(i) = Clw(i) + dt*(-Ux*(Clw(i+1)-Clw(i-1))/(2*dx) +...
               Ux^2*dt/2*(Clw(i+1)-2*Clw(i)+Clw(i-1))/dx^2 +...
               D*(Clw(i+1)-2*Clw(i)+Clw(i-1))/dx^2); % Lax-Wendroff
 end
 Clwnxp1 = Clw(2);
 i=Nx;
 Clwnew(i) = Clw(i) + dt*(-Ux*(Clwnxp1-Clw(i-1))/(2*dx) +...
             Ux^2*dt/2*(Clwnxp1-2*Clw(i)+Clw(i-1))/dx^2 +...
             D*(Clwnxp1-2*Clw(i)+Clw(i-1))/dx^2); % Lax-Wendroff 
  

Same for CU, I guess.

[Logan Parrish]

Thanks for the help, it cleared up a lot of things.