2
$\begingroup$

I have $v_t+v_x=0$ with initial condition $v(x,0)= \sin^2 \pi(x-1)$ for $ x \in [1,2]$. My goal is to find numerical solution for $x \in [0,8]$ using the Leapfrog scheme

$$ u_k^{n+1} = u_k^{n-1} - \frac{ \Delta t }{\Delta x} (u_{k+1}^n - u_{k-1}^n ) $$

I implemented this on matlab using $100$ nodal points and $\Delta t = .75 \Delta x $. Here is my code:

clc;
clear;

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

xmin=0;
xmax=8;
N=100;
dx=(xmax-xmin)/N; %%number of nodes-1

t=0;
tmax=4;
dt=0.75*dx;
tsteps=tmax/dt;

%%%%%Discretization%%%%%%%%
x=xmin-dx:dx:xmax+dx;

%%%%%initial conditions
u0 = F(x);

%%%%Here I am finding u_k^1, I used Euler's method to do so.
u1 = zeros(1,11);
    for k=2:N+1
        u1(k)=u0(k)-(dt/dx)*(u0(k+1)-u0(k));
        u1(1)=0;
        u1(N+2)=0;
        u1(N+3)=0;
    end

u=u0;
unew = u1;

for n=1:tsteps

    u(1)=u(3);
    u(N+3)=u(N+1);

    for i=2:N+2
        unew(i)=u0(i)-(dt/dx)*(u(i+1)-u(i-1));
    end

    t=t+dt;
    u=unew;
    plot(x,u,'bo-')

end

Here is the plot of $u$ vs $x$ when $t=4$.

enter image description here

Which is blown up. Is my code correct? Any criticism would be greatly appreciated.

$\endgroup$
3
$\begingroup$

The implementation of two-step methods such as the Leapfrog scheme requires more data vectors than the implementation of one-step methods. One should not forget to update all the data vectors while iterating. Lastly, the initialization must be performed in an upwind fashion. Here is a revised code:

%%%%Here I am finding u_k^1, I used Euler's method to do so.
u1 = u0;
for k=2:N+1
    u1(k)=u0(k)-(dt/dx)*(u0(k)-u0(k-1));
end
u1(1)=u1(2);

unm1=u0;
un=u1;
unp1=u1;

for n=1:tsteps
    for i=2:N
        unp1(i)=unm1(i)-(dt/dx)*(un(i+1)-un(i-1));
    end

    t=t+dt;
    unm1=un;
    un=unp1;
    plot(x,un,'bo-')

end

with the following output

output

$\endgroup$
  • $\begingroup$ Is there a way in matlab to find the order of convergence of the scheme? $\endgroup$ – Mikey Spivak Feb 11 at 11:58
  • $\begingroup$ @Neymar It can be done numerically by measuring the error between the analytical solution and the numerical solution for various mesh sizes. $\endgroup$ – Harry49 Feb 11 at 12:30

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.