我如何修复我的代码以找到收敛率

Question

IBVP是

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)

我必须首先使用SPACE中的四阶中央差SBP运算符和RK4及时使用Spacing XI = I*H，i = 0,1，.. 。，n和在每个RK阶段更新g（t）。同样，使用代码计算一系列网格的收敛速率。

下面我展示了我的工作，这并不能帮助我找到融合率，所以任何人都可以帮助我找到错误。

%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**

Answer 1

一个出色的错误是

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

您需要括号，否则您可以计算h为xn-（x0/nstp）-1 = 1-（0/nstp）-1 = 0。

或避免所有这些并使用实际的网格步骤。另外，如果nstp是细分中的步骤或段数，则节点的数量更大。而且，如果您定义x0，xn，则实际上也应该使用它们。

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

---

至于SBP-SAT方法，空间离散化给出了

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

摩尔ode函数，因为必须

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

调整矩阵构造以对这些更改的操作具有正确的扩展。

总共提供了实施RK4，

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

---

在下面给出了一个修改的脚本。稳定的解决方案仅适用于cfl = 1或以下。因此，错误的行为似乎是第四阶。初始条件的因素在0.5至1的范围内也最有效，较大的因素倾向于在x = 0时增加误差，至少在初始时间步骤中。在较小的因素的情况下，误差在整个间隔内是均匀的。

%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**