Lax格式激波管程序

clear
xmax=10;
dx=xmax/10000;
x=0:dx:xmax;
M=length（x）;
timemax=0.0045;
dt=1e-6;
gama=1.4;R=287;
for i=1:M
    T（i）=300;
    u0（i）=0;
end
for i=1:5000
    p0（i）=1e6;
    rho0（i）=p0（i）/（R*T（i））;
    E0（i）=R*T（i）/（gama-1）+0.5*u0（i）^2;
    H0（i）=gama*R*T（i）/（gama-1）+0.5*u0（i）^2;
end
for i=5001:M
    p0（i）=1e5;
    rho0（i）=p0（i）/（R*T（i））;
    E0（i）=R*T（i）/（gama-1）+0.5*u0（i）^2;
    H0（i）=gama*R*T（i）/（gama-1）+0.5*u0（i）^2;
end
time=0;
while time    time=time+dt;
    for i=2:M-1
        rho（i）=0.5*（rho0（i+1）+rho0（i-1））-（rho0（i+1）*u0（i+1）-rho0（i-1）*u0（i-1））*dt/（2*dx）;%连续性方程
        U（i）=0.5*（rho0（i+1）*u0（i+1）+rho0（i-1）*u0（i-1））-（rho0（i+1）*u0（i+1）^2+p0（i+1）-rho0（i-1）*u0（i