WENO5-IMEX
Posted yuewen-chen
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了WENO5-IMEX相关的知识，希望对你有一定的参考价值。
- WENO5-IMEX
- 内层使用5点差分，Newmann boundary condition 使用Hermit interpolation ；
- 外层使用WENO5 , ghost point 使用 Lagrange interpolation； 不施加边界条件；
- IMEX 的 (Delta) 矩阵，一定要使用边界条件，不能使用周期边界条件；
- 这次还用有一个创新，使用了MATHEMATICA 计算插值系数，带来很大方便。
Module[{},
x2 = x1 + h;
x3 = x1 + 2*h;
x4 = x1 + 3*h;
x5 = x1 + 4*h;
eq1[x_] := a[1] + a[2]*x + a[3]*x^2 + a[4]*x^3 + a[5]*x^4;
Th = First@
Solve[eq1[x1] == y[1] && eq1[x2] == y[2] && eq1[x3] == y[3] &&
eq1[x4] == y[4] && eq1[x5] == y[5], {a[1], a[2], a[3], a[4],
a[5]}];
neq[x_] := eq1[x] /. Th // Simplify;]
(neq''[x3] /. {y -> u}) // Simplify;
Module[{},(*Hermit 插值*)
x2 = x1 + h;
x3 = x1 + 2*h;
x4 = x1 + 3*h;
x5 = x1 + 4*h;
eq2[x_] := b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3 + b[5]*x^4;
eq3[x_] := b[2] + 2*b[3]*x + 3*b[4]*x^2 + 4*b[5]*x^3;
Th2 = First@
Solve[eq2[x1] == v[1] && eq2[x2] == v[2] && eq2[x3] == v[3] &&
eq2[x4] == v[4] && eq3[x1] == dv[1], {b[1], b[2], b[3], b[4],
b[5]}];
eq4[x_] := eq2[x] /. Th2 // Simplify;
]
Module[{},
f1 = -((u[2] - 16 u[3] + 30 u[4] - 16 u[p1] + u[p2])/(12 h^2));
f2 = f1 /. {u[p1] -> 1/3 (u[1] - 6 u[2] + 18 u[3] - 10 u[4]),
u[p2] -> (8 u[1])/3 - 15 u[2] + 40 u[3] - (80 u[4])/3} //
Simplify;(*f2 是右端的插值*)
f2 = f2 /. {u[4] -> u[n], u[3] -> u[n - 1], u[2] -> u[n - 2],
u[1] -> u[n - 3]};
f3 = -((u[m2] - 16 u[m1] + 30 u[1] - 16 u[2] + u[3])/(12 h^2));
f4 = f3 /. {u[m1] -> 1/3 (-10 u[1] + 18 u[2] - 6 u[3] + u[4]),
u[m2] -> -((80 u[1])/3) + 40 u[2] - 15 u[3] + (8 u[4])/3} //
Simplify;(*f1 是左端插值*)
];
Clear[u, f5, f6]
Module[{},(*矩阵第二行的插值计算*)
f5 = -((u[m1] - 16 u[1] + 30 u[2] - 16 u[3] + u[4])/(
12 h^2));(*f5 是左端 多出一个的 ghost point*)
f6 = f5 /. {u[m1] -> 1/3 (-10 u[1] + 18 u[2] - 6 u[3] + u[4])} //
Simplify;
f7 = -((u[1] - 16 u[2] + 30 u[3] - 16 u[4] + u[p1])/(
12 h^2));(*f7 是右端 多出一个的 ghost point*)
f8 = f7 /. {u[p1] -> 1/3 (u[1] - 6 u[2] + 18 u[3] - 10 u[4])} //
Simplify;
f8 = f8 /. {u[4] -> u[n], u[3] -> u[n - 1], u[2] -> u[n - 2],
u[1] -> u[n - 3]};(*直接打印出来 f8*)
]
Coefficient[f8, {u[n - 3], u[n - 2], u[n - 1], u[n]}]function T80
N=200;
eps=1/20;
a=5;
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:2)),df(x(end-1:end))];
%u=[An(1:30);An(31:70);An(71:100)]*0.8;
u=sin(x);
t=0;
dt=0.5*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=D1_5points(x,u,N_BC,'N');
p=du./sqrt(1+d1u.^2);
%-----2. 外部使用WENO5---------
d2u_p=WENO5_1D(x,p,D_BC, 1,'no','no');
d2u_m=WENO5_1D(x,p,D_BC,-1,'no','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]);
um2=circshift(u,[2,0]);
up1=circshift(u,[-1,0]);
up2=circshift(u,[-2,0]);
g1=-(80*u(1))/3.+40*u(2)-15*u(3)+(8*u(4))/3.;
g2=(-10*u(1) + 18*u(2) - 6*u(3) + u(4))/3.;
g3=(u(end-3)-6*u(end-2)+18*u(end-1)-10*u(end))/3.;
g4=(8*u(end-3))/3.-15*u(end-2)+40*u(end-1)-(80*u(end))/3.;
um1(1)=g2;
um2(1)=g1;
um2(2)=g2;
up1(end) =g3;
up2(end-1)=g3;
up2(end) =g4;
y=-1/(12*h^2)*(um2-16*um1+30*u-16*up1+up2);
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=-16*ones(n-1,1);
c3= 30*ones(n, 1);
M=sparse(diag(c1,-2)+diag(c2,-1)+diag(c3,0)+diag(c2,1)+diag(c1,2))/(-12*h^2);
M(1,1)=-85/(18*h^2);
M(1,2)=108/(18*h^2);
M(1,3)=-27/(18*h^2);
M(1,4)=4/(18*h^2);
M(n,n-3)= 4/(18*h^2);
M(n,n-2)=-27/(18*h^2);
M(n,n-1)=108/(18*h^2);
M(n,n) =-85/(18*h^2);
M(2,1:4)=[29/(18*h^2),-(3/(h^2)), 3/(2*h^2), -(1/(9*h.^2))];
M(n-1,n-3:n)=[-(1/(9*h^2)), 3/(2*h^2), -(3/h^2), 29/(18*h^2)];
y=M;
end
end以上是关于WENO5-IMEX的主要内容，如果未能解决你的问题，请参考以下文章