Hamilton Jacobi

function T77
N=500;
x=linspace(0,4*pi,N)';
h=x(2)-x(1);
N_BC=[cos(x(1:2)),cos(x(end-1:end))];
D_BC=[sin(x(1:2)-2*h),sin(x(end-1:end)+2*h)];
u_0=sin(x);
u=[u_0(1:100);cos(x(101:end-101)).*u_0(101:end-101);u_0(end-100:end)];
t=0;
dt=0.8*h;
t_end=50;
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;
    subplot(1,2,1)
    plot(x,u,'b.',x,sin(x),'r')
    subplot(1,2,2)
    plot(x,err(u),'.')
    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

Hamilton Jacobi

标签：splay   ace   point   cup   sub   play   UNC   lin   差分   

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