二维euler方程组harten数值格式 matlab代码

以下是二维Euler方程组的Harten数值格式的Matlab代码：

% 定义初始条件和网格参数
nx = 100; ny = 100;
x = linspace(0,1,nx); y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,y);
dx = x(2) - x(1); dy = y(2) - y(1);
dt = 0.001; tmax = 0.1; nt = ceil(tmax/dt);

% 定义物理参数
gamma = 1.4;
xmin = min(min(xx)); xmax = max(max(xx));
ymin = min(min(yy)); ymax = max(max(yy));
xc = 0.5*(xmin + xmax); yc = 0.5*(ymin + ymax);
r = sqrt((xx - xc).^2 + (yy - yc).^2);
rho = 1 - 0.3146*r;
p = 1 - 0.1*r;
u = 0.92*ones(size(xx));
v = zeros(size(xx));
E = p./(gamma - 1) + 0.5*rho.*(u.^2 + v.^2);

% 定义数值格式中的常数
beta = 1.5;
eps = 0.1*dx;

% 进行时间循环
for n = 1:nt
    % 计算数值通量
    fx = rho.*u; fy = rho.*v;
    Etotal = E + 0.5*rho.*(u.^2 + v.^2);
    H = Etotal + p;
    c = sqrt(gamma*p./rho);
    lambda1 = u + c; lambda2 = u - c;
    lambda3 = v + c; lambda4 = v - c;
    alpha1 = max(0, -beta*min(lambda1,0), beta*max(lambda1,0));
    alpha2 = max(0, -beta*min(lambda2,0), beta*max(lambda2,0));
    alpha3 = max(0, -beta*min(lambda3,0), beta*max(lambda3,0));
    alpha4 = max(0, -beta*min(lambda4,0), beta*max(lambda4,0));
    f1 = fx; f2 = fx.*u + p; f3 = fy; f4 = fy.*v + p;
    f1bar = 0.5*(f1 + circshift(f1,[0,-1])); f2bar = 0.5*(f2 + circshift(f2,[0,-1]));
    f3bar = 0.5*(f3 + circshift(f3,[-1,0])); f4bar = 0.5*(f4 + circshift(f4,[-1,0]));
    f1tilde = f1bar - 0.5*dt/dx*(alpha1.*circshift(f1 - f1,[0,-1]) - alpha2.*(f1 - circshift(f1,[0,-1])));
    f2tilde = f2bar - 0.5*dt/dx*(alpha1.*circshift(f2 - f2,[0,-1]) - alpha2.*(f2 - circshift(f2,[0,-1]))) ...
                      - dt/dx*(alpha1.*circshift(H - H,[0,-1]) - alpha2.*(H - circshift(H,[0,-1])));
    f3tilde = f3bar - 0.5*dt/dy*(alpha3.*circshift(f3 - f3,[-1,0]) - alpha4.*(f3 - circshift(f3,[-1,0])));
    f4tilde = f4bar - 0.5*dt/dy*(alpha3.*circshift(f4 - f4,[-1,0]) - alpha4.*(f4 - circshift(f4,[-1,0]))) ...
                      - dt/dy*(alpha3.*circshift(H - H,[-1,0]) - alpha4.*(H - circshift(H,[-1,0])));

    % 计算数值解
    rho = rho - dt/dx*(f1tilde - circshift(f1tilde,[0,1])) - dt/dy*(f3tilde - circshift(f3tilde,[1,0]));
    rho = max(rho,eps);
    u = (f2tilde - circshift(f2tilde,[0,1]))./rho;
    v = (f4tilde - circshift(f4tilde,[1,0]))./rho;
    E = Etotal - 0.5*rho.*(u.^2 + v.^2);
    p = (gamma - 1)*rho.*(E - 0.5*(u.^2 + v.^2));
end

% 画图
figure
pcolor(xx,yy,rho); shading interp; colorbar; title('Density')
figure
pcolor(xx,yy,p); shading interp; colorbar; title('Pressure')
figure
quiver(xx,yy,u,v); title('Velocity')