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.
L=2;
D = 0.001;
Nx=50;
Nt=1000;
dt=0.01;
x=linspace(0,L,Nx);
dx=x(2)-x(1);
t=0;
Ci = @(x) cos(2*pi*x/L);
CU = Ci(x);
Clw = Ci(x);
Ux = 0.5;
CUnew = CU;
Clwnew = Clw;
for n=1:Nt
t=t+dt;
for i=2:Nx-1
if Ux>0
dCUdx=(CU(i)-CU(i-1))/dx;
d2CUdx2 = (CU(i-1)-2*CU(i)+CU(i+1))/(dx^2);
else
dCUdx= (CU(i+1)-CU(i))/dx;
d2CUdx2 = (CU(i-1)-2*CU(i)+CU(i+1))/(dx^2);
end
CUnew(i)=CU(i)+dt*(-Ux*dCUdx + D*d2CUdx2);
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));
C_analytic = @(x,t) cos(2*pi*(x(i)-Ux*t(n))/L)*exp(-D(2*pi/L)^2*t(n));
end
CUnew(1) = CUnew(Nx-1);
CUnew(Nx) = CUnew(2);
Clwnew(1) = Clwnew(Nx-1);
Clwnew(Nx) = Clwnew(2);
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
