三阶CWENO的导数逼近求解线性对流方程的Matlab代码

以下是使用三阶CWENO格式求解线性对流方程的Matlab代码：

% 设置参数
c = 1; % 对流速度
dx = 0.1; % 空间步长
dt = 0.01; % 时间步长
x = 0:dx:2; % 空间网格
t = 0:dt:2; % 时间网格
nx = length(x); % 空间网格数
nt = length(t); % 时间网格数
f = @(x) exp(-100*(x-1).^2); % 初始条件

% 初始化解
u = zeros(nx,nt);
u(:,1) = f(x);

% 计算系数
beta0 = 1/3; beta1 = 2/3;

% 迭代求解
for n = 1:nt-1
    % 计算数值通量
    uL = [u(nx,n) u(:,n) u(:,n)];
    uR = [u(:,n) u(:,n) u(1,n)];
    fL = c*(1-beta0)*(uL(:,2:end)-uL(:,1:end-1))/dx + c*beta0*(uL(:,3:end)-uL(:,2:end-1))/dx;
    fR = c*(1-beta0)*(uR(:,2:end)-uR(:,1:end-1))/dx + c*beta0*(uR(:,3:end)-uR(:,2:end-1))/dx;
    fBar = beta1*fL + (1-beta1)*fR;
    % 计算斜率
    ux = zeros(nx+2,1);
    ux(2:end-1) = (u(:,n)-u(:,n-1))/dt;
    uxL = (1-beta0)*(ux(2:end-1)-ux(1:end-2))/dx + beta0*(ux(3:end)-ux(2:end-1))/dx;
    uxR = (1-beta0)*(ux(3:end)-ux(2:end-1))/dx + beta0*(ux(4:end)-ux(3:end))/dx;
    sL = sign(uxL); sR = sign(uxR);
    alphaL = sL.*max(0,abs(uxL)); alphaR = sR.*max(0,abs(uxR));
    wL = alphaL./(alphaL+alphaR); wR = alphaR./(alphaL+alphaR);
    % 计算逼近值
    uBar = zeros(nx+2,1);
    uBar(2:end-1) = wL.*(u(2:end-1,n)+uxL*dx/2) + wR.*(u(3:end,n)-uxR*dx/2);
    % 更新解
    u(:,n+1) = uBar(2:end-1);
end

% 绘制解的图像
[X,Y] = meshgrid(t,x);
surf(X,Y,u');
xlabel('t'); ylabel('x'); zlabel('u');

其中，数值通量使用了Lax-Friedrichs格式，斜率使用了MUSCL格式，逼近值使用了三阶CWENO格式。