Implementing Lax-Wendroff scheme for advection in matlab

Question

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.

Answer 1

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