二维Burgers方程的微分求积法求解（MATLAB代码）

% DQM法计算齐次第一类边界条件的二维Burgers方程
clear all;close all;
format short e
global D1_y D1_x D2_y D2_x v
v=0.01;
% T1=[0.05,0.25,0.5];
T1=[0.2,0.6,1,1.4];
% T1=[0.6];
for k1=1:length(T1)
     T=T1(k1);
% T=0.05;
%%%x方向离散
    a=0;b=1;
%a=-0.5;b=0.5;
    Nx = 11;
 x=cglgrid(a,b,Nx);%（取CGL点）
 [D1_x,D2_x]=weigcoef(x);%权系数矩阵的计算
%%y方向离散
    c=0;d=1;
%  c=-0.5;d=0.5;
    Ny = 50;
 y=cglgrid(c,d,Ny);[D1_y,D2_y]=weigcoef(y);
%%%%%%%%%%时间离散%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dt1=[0.01,0.001,0.0001];
% for k1=1:length(dt1)
%     dt=dt1(k1);
    dt=0.01;
    t=0:dt:T;
    
 %%%%%%%%%%解析解%%%%%%%%%%%%%%   
    [X,Y]=meshgrid(x,y);
    for k=1:size(t,2)
          u(:,:,k)=1./(1+exp((X+Y-t(k))/(2*v)));
%         u(:,:,k)=0.5-tanh((X+Y-t(k))/(2*v));
    end
    %%%%%%%%%%%%%%%%%%%%
    % 图示解析解
 subplot(121)
surf(x,y,u(:,:,end));
% %     contourf(x,y,u(:,:,end));
 set(gca,'YDir','reverse');
 colorbar;
 title(['解析解: t=', num2str((length(t)-1)*dt)]);
shading flat
    xlabel('x'); ylabel('y');zlabel('u');
    %%%%%%初始条件的输入%%%%%%%%%%%
    uu(:,:,1)=u(:,:,1);
    %求解
    for k=2:length(t)
        %向前欧拉
%         uu(:,:,k)=uu(:,:,k-1)-dt*uu(:,:,k-1).*(D1_y*uu(:,:,k-1))...
%             -dt*uu(:,:,k-1).*(uu(:,:,k-1)*(D1_x)')+dt*v*(D2_y*uu(:,:,k-1)+...
%             uu(:,:,k-1)*(D2_x)');
        %SSPRK43
        uu1=uu(:,:,k-1)+0.5*dt*L(uu(:,:,k-1));
        uu2=uu1+0.5*dt*L(uu1);
        uu3=(2/3)*uu(:,:,k-1)+(1/3)*uu2+(1/6)*dt*L(uu2);
        uu(:,:,k)=uu3+0.5*dt*L(uu3);
%         %%%加边界条件
        uu(1,:,k)=u(1,:,k);uu(end,:,k)=u(end,:,k);
        uu(:,1,k)=u(:,1,k);uu(:,end,k)=u(:,end,k);
    end
%%%%%%%%%%%%%
subplot(122)
    surf(x,y,uu(:,:,end));
    set(gca,'YDir','reverse');
    colorbar;shading flat
    title(['近似解: t=',num2str((k-1)*dt)]);
    xlabel('x');ylabel('y');zlabel('u');
    error=abs(uu(:,:,end)-u(:,:,end));
    error_infty(k1,1)=max(max(error))
    error_L2(k1,1)=sqrt(sum(sum(error.^2))/(Ny*Nx))
end