WENO5-IMEX

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

标签：code   矩阵   一个   length   point   int   new   alc   plot   

原文地址：https://www.cnblogs.com/yuewen-chen/p/11641811.html