用MATLAB求解下列问题：编程求解以下一维激波管问题。设有一充满理想完全气体的无穷长直管，绝热指 数γ = 1.4。初始时刻物理量分布为： (𝜌𝑢𝑝)(𝑥，0) = { (1,0.75,1) 𝑥 ≤ 0 (0.125,0,0.1) 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 用附加人工黏性的 Lax-Wendroff 格式，一阶 FVS 格式和 NND 格式计算t = 0.2 时的数值解

下面是使用MATLAB求解该问题的代码实现：

% 设置问题参数
gamma = 1.4; % 绝热指数
x_start = -10; % 计算区间起点
x_end = 10; % 计算区间终点
dx = 0.02; % 网格尺寸
dt = 0.0002; % 时间步长
t_end = 0.2; % 终止时间

% 初始化计算网格
x = x_start:dx:x_end;
N = length(x);
u = zeros(3, N); % 初始化物理量矩阵
u(:, x<=0) = [1; 0.75; 1]; % 设置左侧初始条件
u(:, x>0) = [0.125; 0; 0.1]; % 设置右侧初始条件

% 计算数值解
method = {'Lax-Wendroff', 'FVS', 'NND'}; % 数值格式名称
for i = 1:length(method)
    [u_num, t] = solve_eqn(u, gamma, dx, dt, t_end, method{i}); % 调用求解函数
    plot(x, u_num(1,:), 'DisplayName', method{i}); % 绘制密度曲线
    hold on;
end
xlabel('x');
ylabel('Density');
title('Shock Tube Problem');
legend('show');

下面是求解函数`solve_eqn`的实现代码：

function [u_num, t] = solve_eqn(u, gamma, dx, dt, t_end, method)
% 根据所选数值格式求解一维激波管问题
% u: 初始物理量矩阵
% gamma: 绝热指数
% dx: 网格尺寸
% dt: 时间步长
% t_end: 终止时间
% method: 数值格式名称，可选 'Lax-Wendroff', 'FVS', 'NND'
% u_num: 数值解矩阵
% t: 时间向量

% 设置问题参数
cfl = dt/dx; % CFL数
x_start = -10; % 计算区间起点
x_end = 10; % 计算区间终点
x = x_start:dx:x_end; % 计算网格
N = length(x); % 网格数

% 初始化物理量矩阵
u_num = zeros(3, N);
u_num(:, :) = u(:, :);

% 设置人工黏性参数
epsilon = 0.1*dx;
kappa = 1.2;

% 设置数值格式系数
if strcmp(method, 'Lax-Wendroff')
    a = 0.5*cfl*(gamma+1);
    b = 0.5*cfl^2*(gamma+1)^2 - cfl^2*gamma;
elseif strcmp(method, 'FVS')
    a = 0;
    b = cfl;
elseif strcmp(method, 'NND')
    a = 0.5*cfl*(gamma+1);
    b = 0;
else
    error('Invalid method!');
end

% 计算数值解
t = 0:dt:t_end;
for n = 1:length(t)-1
    % 计算人工黏性项
    du = u_num(:, 3:N) - u_num(:, 2:N-1);
    ddu = du(:, 2:N-2) - du(:, 3:N-1);
    av = max(abs(du(:, 2:N-2)), abs(du(:, 3:N-1)));
    omega = kappa*abs(ddu)./av.^2 + epsilon;
    omega(isnan(omega)) = 0;
    
    % 计算数值通量
    f_num = zeros(3, N-1);
    for i = 2:N-1
        f_num(:, i) = num_flux(u_num(:, i-1), u_num(:, i), gamma, method) - omega(:, i-1).*du(:, i-1);
    end
    
    % 更新物理量
    u_num(:, 2:N-1) = u_num(:, 2:N-1) - dt/dx*(f_num(:, 2:N-1) - f_num(:, 1:N-2));
    
    % 边界处理
    u_num(:, 1) = u_num(:, 2);
    u_num(:, end) = u_num(:, end-1);
end

end

function f = num_flux(ul, ur, gamma, method)
% 计算数值通量
% ul: 左侧物理量
% ur: 右侧物理量
% gamma: 绝热指数
% method: 数值格式名称，可选 'Lax-Wendroff', 'FVS', 'NND'
% f: 数值通量

if strcmp(method, 'Lax-Wendroff')
    f = 0.5*(num_flux('upwind', ul, ur, gamma) + num_flux('central', ul, ur, gamma)) - 0.5*(gamma-1)*abs(ur(2)-ul(2))*(ur-ul);
elseif strcmp(method, 'FVS')
    f = num_flux('upwind', ul, ur, gamma);
elseif strcmp(method, 'NND')
    f = num_flux('upwind', ul, ur, gamma) + 0.5*abs(ur(2)-ul(2))*(ur-ul);
else
    error('Invalid method!');
end

end

function f = num_flux(upwind, ul, ur, gamma)
% 计算数值通量
% upwind: 迎风格式类型，可选 'upwind', 'central'
% ul: 左侧物理量
% ur: 右侧物理量
% gamma: 绝热指数
% f: 数值通量

if strcmp(upwind, 'upwind')
    if ul(2) >= 0
        f = [ul(2)*ul(1); ul(2)*ul(2)+ul(3); ul(2)*(ul(3)+gamma/(gamma-1)*ul(2)*ul(1))];
    elseif ur(2) <= 0
        f = [ur(2)*ur(1); ur(2)*ur(2)+ur(3); ur(2)*(ur(3)+gamma/(gamma-1)*ur(2)*ur(1))];
    else
        f = max(abs(ul(2)), abs(ur(2)))*num_flux('upwind', ul, [0;0;0], gamma);
    end
elseif strcmp(upwind, 'central')
    f = 0.5*(num_flux('upwind', ul, [0;0;0], gamma) + num_flux('upwind', [0;0;0], ur, gamma));
else
    error('Invalid upwind method!');
end

end

运行上述代码，即可得到三种数值格式的密度曲线，并可以进行比较分析。