Test code for 1D

  • Test code for 1D

    In order to test the WENO5 method and the 5 point center difference method, we test follow simple case

\[\begin{align} u_t&=u_{xx}+\sin(x),x\in [0,2\pi]\\ u_x|_{x=0} &=\sin(0)\\ u_x|_{x=2\pi} &=\sin(2\pi)\\ \end{align} \]

we treat \(u_x\) with 5 points center difference , no boundary condition, the outside use Dirichlet boundary condition.

  • Step1 . we use implicit scheme, RungeKutta-WENO5 method.
    The optimal third order TVD Runge-Kutta method is given by follow:

\[ \begin{align} u_t &=L(u)\\ u_1 &=u^n+dt L(u^n)\\ u_2 &=\frac{3}{4}u^n+\frac{1}{4}u_1+\frac{1}{4}dtL(u_1)\\ u^{n+1} &=\frac{1}{3}u^n+\frac{2}{3}u_2+\frac{2}{3}dtL(u_2) \end{align} \]

发现了一些而问题:WENO5 的Dirichlet 边界条件施加的时候,不需要插值,直接放置左右各两点即可, 启用 ‘smooth’ 模式比较好。内部五点中心差分不需要边界条件,对ghost points 使用 Lagrange5 插值即可。
现在任然要调查, Jang equation 出现的边界条件不符合的问题,问题可能出现在三个部分:1. 二阶部分:内部中心差分,外部WENO5, 刚刚排除了这种可能,使用内部无边界条件,外部‘smooth’模式,是成功的,对上面的Possion equation 是好的;2. Hamilton-Jacobi 部分,这需要检验一下。3.第三部分,就是IMEX 的部分,使用了循环边界条件,但是根据以前的经验,这个似乎是没有影响的。

Step2. test Hamilton-Jacobi equation part

\[ \begin{align} u_t-\frac{u_x^2}{1+u_x^2} &=-\frac{\cos(x)^2}{1+\cos(x)^2}\\ u_x&=NewBC \end{align} \]

有问题的,以前的边界条件提的太简单,问题没有出现,现在发现问题不这么简单。\
为了测试是不是 WENO5 造成的,改成差分:

\[\begin{align} u_i^c&=\frac{1}{12 h}(u_{i-2}-8 u_{i-1}+8 u_{i+1}-u_{i+2})\\ u_{i}^+ &=-\frac{1}{12 h}(u_{i-3}-6 u_{i-2}+18 u_{i-1}-10 u_i-3 u_{i+1})\\ u_{i}^- &=\frac{1}{12 h}(-3 u_{i-1}-10 u_i+18 u_{i+1}-6 u_{i+2}+u_{i+3}) \end{align} \]

使用了差分法任然和WENO5 是一样的结果,那么身下的就是边界的问题或者是 flux 的问题,使用 Lax Fridrich 格式,还是一样,所以排除 flux 的问题,剩下的就是边界问题了。
可以考虑inverse Lax Wendroff 格式了。

function T76
N=100;
x=linspace(0,2*pi,N)';
h=x(2)-x(1);
D_BC=[cos(x(1)-2*h),cos(x(1)-h),cos(x(end)+h),cos(x(end)+2*h)];
N_BC=[sin(x(1)),sin(x(end))];
U1=sin(x);
U2=x.^2;
%u=[U1(1:5);U2(6:end-11);U1(end-10:end)];
u=cos(x)*0.1;
t=0;
dt=h^2;
t_end=20;
%============= Runge-Kutta =================
while t<t_end
    u1=u+dt*L(u);
    u2=3/4*u+1/4*u1+1/4*dt*L(u1);
    u=1/3*u+2/3*u2+2/3*dt*L(u2);
    t=t+dt;
     plot(x,u,x,sin(x),'r')
    title(['t=',num2str(t)])
    drawnow
end

   function y=L(u)
        %-----1.内部使用中心差分-------
         du=D1_5points(x,u,N_BC,'no');
        %-----2. 外部使用WENO5---------
         d2u_p=WENO5_1D(x,du,D_BC, 1,'D','smooth');
         d2u_m=WENO5_1D(x,du,D_BC,-1,'D','smooth');
         d2u=(d2u_p+d2u_m)/2;
         y=d2u+sin(x);
    end
end
function T77
N=40;
x=linspace(pi/2,3/2*pi,N)';
h=x(2)-x(1);
N_BC=[cos(x(1:2)),cos(x(end-1:end))];
u=sin(x)*0.8;
t=0;
dt=0.5*h^2;
t_end=20;
S=cos(x).^2./(1+cos(x).^2);
%============= Runge-Kutta =================
while t<t_end
    u1=u+dt*L(u);
    u2=3/4*u+1/4*u1+1/4*dt*L(u1);
    u=1/3*u+2/3*u2+2/3*dt*L(u2);
    t=t+dt;
    
    plot(x,u,'b.',x,sin(x),'r')
    title(['t=',num2str(t)])
    drawnow
end


function y=L(u)
      dup=WENO5_1D(x,u,N_BC, 1,'N','no');
      dum=WENO5_1D(x,u,N_BC,-1,'N','no');
     %-------- 这是差分法 ------------
       %dup=FD5_point(x,u,N_BC,-1,'N'); 
       %dum=FD5_point(x,u,N_BC,1,'N'); 
      
      %y=-Ham(dup,dum)+S;  
       y=-Lax_Fridrich(dup,dum)+S;
       
end

function y=Ham(fxp,fxm)
               by=(min(fxp,0)).^2+(max(fxm,0)).^2;
               y=by./(1+by);
end

    function y=Lax_Fridrich(dup,dum)
            du=(dup+dum)/2;
          alpha=max((-2*power(du,3))./power(1 + power(du,2),2) + (2*du)./(1 + power(du,2)));
          y=du.^2./(1+du.^2)-alpha*(dup-dum);
    end

    function y=err(u)
        du=WENO5_1D(x,u,N_BC, 1,'N','smooth');
        Rh=du.^2./(1+du.^2);
        y=Rh-S;
    end
end
posted @ 2019-10-02 19:14  yuewen_chen  阅读(179)  评论(0编辑  收藏  举报