0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ What is you error ? $\endgroup$ – Paul Cottalorda Feb 4 at 0:43
  • $\begingroup$ 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); $\endgroup$ – Mikey Spivak Feb 4 at 0:57
  • $\begingroup$ 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. $\endgroup$ – Paul Cottalorda Feb 4 at 1:04
  • $\begingroup$ 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\}$. $\endgroup$ – Christoph Feb 4 at 3:28
  • $\begingroup$ Additionally, your value of $p=2$ is too large, you need $p\le 1$ to get a correct propagation of the wave. $\endgroup$ – LutzL Feb 4 at 15:12
4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.