Implementing Lax-Wendroff scheme for advection in matlab

MATLABnumerical methodspartial differential equationstransport-equation

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.

Best 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 enter image description here

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
Related Question