- WENO3-IMEX
- 内层使用中心差分,边界点使用 Hermit 插值,施加 Newmann 边界条件;
- 外层使用WENO3, 不施加边界条件,ghost point 使用WENO3 外插值;
- IMEX 操作的时候,都要施加边界条件。前几次的失败都是相信,使用循环边界条件也可以,实际上不行!
- 假如ghost point 使用 lagrange 插值会怎样?
function T82
N=200;
eps=1/10;
a=2 ;
x=linspace(0,2*pi,N)';
r0=(x(1)+x(end))/2;
x0=r0;
h=x(2)-x(1);
df=@(r)1./(power(eps,2) + power(-r0 + r,2));
mean_df=@(r)df(r)./sqrt(df(r).^2+1);
d1u=df(x);
ss=@(x)((-2*x + 2*x0)./(sqrt(1 + power(power(eps,2) + power(x - x0,2),-2)).*...
(1+power(eps,4)+power(x,4) + 2*power(eps,2).*power(x - x0,2) - 4*power(x,3)*x0 +...
6*power(x,2).*power(x0,2) - 4*x.*power(x0,3) + power(x0,4))));
An=-(atan((r0 - x)/eps)/eps);
S1=-ss(x);
D_BC=[mean_df(x(1:2)-2*h),mean_df(x(end-1:end)+2*h)];
N_BC=[df(x(1)),df(x(end))];
u=sin(x);
t=0;
dt=2*h;
t_end=2000;
L=eye(N,N)*(1/dt)-a*D2(x);
%--------------------------------------------
while t<t_end
f_old=L\(S(u));
u=f_old;
t=t+dt;
plot(x,u,'b-.',x,An,'r')
title(['t=',num2str(t)])
drawnow
end
%--------------------------------------------
function s=S(u)
du=D3(x,u,N_BC,'N');
p=du./sqrt(1+d1u.^2);
%-----2. 外部使用WENO5---------
d2u_p=WENO3(x,p,1,'no');
d2u_m=WENO3(x,p,-1,'no');
J=(d2u_p+d2u_m)/2;
s=u/dt+J+S1-a*d2u(x,u);
end
%--------------------------------------------
function y=d2u(x,u)
% The goal of this function is calculate the second order
% derivative of u, with a periodic boundary condition
h=x(2)-x(1);
um1=circshift(u,[1,0]);
up1=circshift(u,[-1,0]);
up1(end)=2*h*N_BC(2)*0 + u(end-1);
um1(1) =-2*h*N_BC(1)*0 + u(2);
y=(um1-2*u+up1)/(h^2);
end
function y=D2(x)
% out put: Matrix D2, with periodic boundary condition
h=x(2)-x(1);
n=length(x);
c1=ones(n-2,1);
c2=ones(n-1,1);
c3=-2*ones(n,1);
M=sparse(diag(c2,-1)+diag(c3,0)+diag(c2,1));
M=M+sparse([1,n],[2,n-1],[1,1],n,n);
y=M/(h^2);
end
end