# Implementing Lax-Wendroff scheme for advection in matlab

Given the advection equation $$v_t + v_x = 0$$ with initial condition $$u(x,0) = \sin^2 \pi(x-1)$$ for $$1 \leq x \leq 2$$ and $$0$$ otherwise. Solve this PDE exactly and compare with numerical solution using the following Lax-Wendroff scheme $$\frac{u_j^{n+1} - u_j^n}{\Delta t} + a \frac{u_{j+1}^n-u_{j-1}^n}{2\Delta x} - \frac{1}{2} a^2 \Delta t \left(\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{\Delta x^2}\right) = 0$$

For the code, take $$0 \leq x \leq 6$$, from $$t=0$$ to $$t=4$$ and $$\Delta t = 0.01$$ and mesh intervals $$N=200$$ (201 grid points)

### My solution

By the method of charactheristic, we see that the solution is $$v(x,t) = F(x-t)$$ since $$v(x,0)= \sin^2 \pi(x-1) = F(x)$$, then $$\boxed{ v(x,t) = \sin^2 \pi(x-t-1)}$$

Now, to implement this scheme in matlab, I rewrite the equation as follows after multiplying by $$\Delta t$$ throughout

$$u_j^{n+1} = u_j^n (1+p^2) + u_{j+1}^n (0.5p - 0.5 p^2) - u_{j-1}^n (0.5(p+p^2))$$

where $$p = \dfrac{ \Delta t }{\Delta x}$$. Here is the code in matlab.

clear;

%%%%Initial conditions, initialization %%%%%%%

F = @(x) sin(pi*(x-1)).^2 .* (1<x).*(x<2);

m=201;
x = linspace(0,6,m);
dx=6/(m-1);
t=0;
dt=0.01;
p=dt/dx;
v=F(x);

figure;
plot(x,v);

%%Here goes the Scheme iteration%%%%
while t<4
v(2:m-1) = (1+p^2)*v(1:m-1)+(p./2-p^2./2)*v(1:m)-(p./2+p^2./2)*v(1:m-2);
v(1)=v(2);
v(m)=v(m-1);
t=t+dt;
end

plot(x,v,'bo');
hold on
plot(x,F(x-t),'k-');


However, I keep getting an error. What is wrong with my code? any help would be greatly appreciated.

• What is you error ? – Paul Cottalorda Feb 4 at 0:43
• Matrix dimensions must agree. Error in LaxWend (line 20) v(2:m-1) = (1+p^2)*v(1:m-1)+(p./2-p^2./2)*v(1:m)-(p./2+p^2./2)*v(1:m-2); – Mikey Spivak Feb 4 at 0:57
• Isn't it because v(1:m), v(1:m-1), v(1:m-2) are not of the same size ? You must add offsets on the left too. – Paul Cottalorda Feb 4 at 1:04
• Indeed, as Paul Cottalorda wrote. If $j \in \{2, \dots, m-1\}$, then $j+1 \in \{3, \dots, m\}$ and $j-1 \in \{1, \dots, m-2\}$. – Christoph Feb 4 at 3:28
• Additionally, your value of $p=2$ is too large, you need $p\le 1$ to get a correct propagation of the wave. – LutzL Feb 4 at 15:12

## 1 Answer

You transformed the method equation wrongly, there are multiple sign errors. It should be \begin{align} u_j^{n+1} &=u_j^n - \frac{p}2(u_{j+1}^n - u_{j-1}^n) + \frac{p^2}{2}(u_{j+1}^n - 2u_j^n + u_{j-1}^n) \\ &= u_j^n (1-p^2) -\frac{p-p^2}2 u_{j+1}^n +\frac{p+p^2}2 u_{j-1}^n \end{align}

Implementing this results in a correct looking dynamic F = @(x) sin(pi*(x-1)).^2 .* (1<x).*(x<2);

n = 200;
x0 = 0; xn = 6;
x = linspace(x0,xn,n+1);
dx = x(2)-x(1);
t = 0;
dt = 0.1*dx;
p = dt/dx;
v = F(x);

%%% Scheme iterations
k=1
while t<4
v(2:n) = (1-p^2)*v(2:n)-0.5*(p-p^2)*v(3:n+1)+0.5*(p+p^2)*v(1:n-1);
v(1) = v(2);
v(n+1) = v(n);
t = t + dt;
if t>k
subplot(4,1,k);
hold on
plot(x,v,'bo');
plot(x,F(x-t),'k-');
ylim([-0.1,1.1]);
grid on;
title(sprintf("t=%.3f",t))
hold off
k=k+1
end
end