How do I fix my code to find the convergence rate

Question

The IBVP is

u_t+2u_x=0 , 0<x<1, t>0

u(x,0)=sin(pi*x)
u(0,t)=g(t)=sin(-2*pi*t)

I have to first implement the scheme using 4th-order central difference SBP operator in space and RK4 in time with spacing xi=i*h, i=0,1,...,N and update g(t) at every RK stage. Also, using the code compute the convergence rate on a sequence of grids.

Below I have shown my working which is not helping me find the convergence rate so can anyone help me with finding my mistakes.

%Parameters
nstp = 16; %number of grid points in space
t0 = 0; %initial time
tend = 1; %end time
x0 = 0; %left boundary
xN = 1; %right boundary
x = linspace(0,1,nstp); %points in space
h = xN-x0/nstp-1; %grid size
cfl = 4;
k = cfl*h; %length of time steps
N = ceil(tend/k); %number of steps in time
k = tend/N; %length of time steps
u0 = sin(pi*x); %initial data
e = zeros(nstp);
e(1) = 1;
e0 = e(:,1);
m = zeros(nstp);
m(1) = sin(-2*pi*tend);
g = m(:,1);

%4th order central SBP operator in space
m=10; %points

H=diag(ones(m,1),0);
H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]);
H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48])));

H=H*h;

HI=inv(H);   
  
D1=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)- ...
    8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));

D1(1:4,1:6)=[-24/17,59/34,-4/17,-3/34,0,0; -1/2,0,1/2,0,0,0; 
4/43,-59/86,0,59/86,-4/43,0; 3/98,0,-59/98,0,32/49,-4/49];
D1(m-3:m,m-5:m)=flipud( fliplr(-D1(1:4,1:6)));

D1=D1/h;

%SBP-SAT scheme
u = -2*D1*x(1:N-1)'-(2*HI*(u0-g)*e);


%Runge Kutta for ODE
for i=1:nstp %calculation loop
    t=(i-1)*k;
    k1=D*u;
    k2=D*(u+k*k1/2);
    k3=D*(u+k*k2/2);
    k4=D*(u+k*k3);
    u=u+(h*(k1+k2+k3+k4))/6; %main equation

    figure(1)
    plot(x(1:N-1),u); %plot
    drawnow
end



%error calculcation
ucorrect = sin(pi*(x-2*tend)); %correct solution
ucomp = u(:,end); %computed solution
errornorm = sqrt((ucomp-ucorrect)'*H*(ucomp-ucorrect)); %norm of error**

I have the impression that the implementation of the scheme is wrong. You are mixing, as far as I can see, parts of a SBP-SAT in time direction into the intended SBP-SAT in space direction. The stage computations should depend on `t`. It often helps to encapsulate the DE part in `du/dt=F(u,t)`, that is, implement `F` as a function, and use only `F` in the RK steps. — Lutz Lehmann, Apr 23 '22 at 09:23

Lutz Lehmann · Answer 1 · 2022-04-23T12:06:54.977

One egregious error is

h = xN-x0/nstp-1; %grid size

You need parentheses here, else you compute h as xN-(x0/nstp)-1=1-(0/nstp)-1=0.

Or avoid this all and use the actual grid step. Also, if nstp is the number of steps or segments in the subdivision, the number of nodes is one larger. And if you define x0,xN then you should also actually use them.

x = linspace(x0,xN,nstp+1); %points in space
h = x(2)-x(1); %grid size

As to the SBP-SAT approach, the space discretization gives

u_t(t) + 2*D_1*u(t) = - H_inv * (dot(e_0,u(t)) - g(t))*e_0

which should give the MoL ODE function as

F = @(u,t) -2*D1*u - HI*e0*(u(1)-g(t))

The matrix construction has to be adapted to have the correct extends for these changed operations.

The implement RK4 as

k1 = F(u,t)
k2 = F(u+0.5*h*k1, t+0.5*h)
...

In total this gives a modified script below. Stable solutions only occur for cfl=1 or below. With that the errors seemingly behave like 4th order. The factor for the initial condition also works best in the range 0.5 to 1, larger factors tend to increase the error at x=0, at least in the initial time steps. With the smaller factors the error is uniform over the interval.


%Parameters
xstp = 16; %number of grid points in space
x0 = 0; %left boundary
xM = 1; %right boundary
x = linspace(x0,xM,xstp+1)'; %points in space
h = x(2)-x(1); %grid size

cfl = 1;

t0 = 0; %initial time
tend = 1; %end time
k = cfl*h; %length of time steps
N = ceil(tend/k); %number of steps in time
k = tend/N; %length of time steps

u = sin(pi*x); %initial data
e0 = zeros(xstp+1,1);
e0(1) = 1;
g = @(t) sin(-2*pi*t);

%4th order central SBP operator in space
M=xstp+1; %points

H=diag(ones(M,1),0); % or eye(M)
H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]);
H(M-3:M,M-3:M)=fliplr(flipud(diag([17/48 59/48 43/48 49/48])));

H=H*h;
HI=inv(H);   
  
D1=(-1/12*diag(ones(M-2,1),2)+8/12*diag(ones(M-1,1),1)- ...
    8/12*diag(ones(M-1,1),-1)+1/12*diag(ones(M-2,1),-2));

D1(1:4,1:6) = [ -24/17,  59/34, -4/17,  -3/34,     0,     0; 
                  -1/2,      0,    1/2,     0,     0,     0; 
                  4/43, -59/86,      0, 59/86, -4/43,     0; 
                  3/98,      0, -59/98,     0, 32/49, -4/49  ];
                  
D1(M-3:M,M-5:M) = flipud( fliplr(-D1(1:4,1:6)));

D1=D1/h;

%SBP-SAT scheme
F = @(u,t) -2*D1*u - 0.5*HI*(u(1)-g(t))*e0;

clf;

%Runge Kutta for ODE
for i=1:N %calculation loop
    t=(i-1)*k;
    k1=F(u,t);
    k2=F(u+0.5*k*k1, t+0.5*k);
    k3=F(u+0.5*k*k2, t+0.5*k);
    k4=F(u+k*k3, t+k);

    u=u+(k/6)*(k1+2*k2+2*k3+k4); %main equation

    figure(1)
    subplot(211);
    plot(x,u,x,sin(pi*(x-2*(t+k)))); %plot
    ylim([-1.1,1.1]);
    grid on;
    legend(["computed";"exact"])
    drawnow
    subplot(212);
    plot(x,u-sin(pi*(x-2*(t+k)))); %plot
    grid on;
    hold on;
    drawnow
end



%error calculcation
ucorrect = sin(pi*(x-2*tend)); %correct solution
ucomp = u; %computed solution
errornorm = sqrt((ucomp-ucorrect)'*H*(ucomp-ucorrect)); %norm of error**

How do I fix my code to find the convergence rate

1 Answers1