一维喷管流动

clc
clear

L=3;               % length of nozzle
nx=31;             % number of grid points
dx=L/（nx-1）;       % mesh size  
C=0.5;             % courant number
shr=1.4;           % specific heat ratio
N=1000;            % number of time steps
x=0:dx:L;          % mesh
A=zeros（1nx）;     % section area
rho=zeros（Nnx）;   % density
T=zeros（Nnx）;     % nondimensional temperature
V=zeros（Nnx）;     % nondimensional velocity
P=zeros（Nnx）;     % nondimensional pressure
Ma=zeros（Nnx）;    % Mach number

% initial conditions
for j=1:nx
    A（j）=1+2.2*（x（j）-1.5）^2;
    rho（1j）=1-0.3146*x（j）;
    T（1j）=1-0.2314*x（j）;
    V（1j）=（0.1+1.09*x（j））*T（1j）^0.5;
    P（1j）=rho（1j）*T（1j）;
    Ma（1j）=V（1j）./（T（1j）.^0.5）;
end


for i=1:N
    % time step
    temp_dt=zeros（1nx）;
    for ii=1:nx
        temp_dt（ii）=C*dx/（V（iii）+T（iii）^0.5）;
    end
    dt=min（temp_dt）;       
    
    % prediction step
    d_rho=zeros（2nx-1）;
    d_V=zeros（2nx-1）;
    d_T=zeros（2nx-1）;
    pre_rho=zeros（nx-1）;
    pre_V=zeros（nx-1）;
    pre_T=zeros（nx-1）;
    for m=1:nx-1
        d_rho（1m）=-V（im）*（rho（im+1）-rho（im））/dx-rho（im）*（V（im+1）-V（im））/dx-rho（im）*V（im）*（log（A（m+1））-log（A（m）））/dx;
        d_V（1m）=-V（im）*（V（im+1）-V（im））/dx-（1/shr）*（（T（im+1）-T（im））/dx+（T（im）/rho（im））*（rho（im+1）-rho（im））/dx）;
        d_T（1m）=-V（im）*（T（im+1）-T（im））/dx-（shr-1）*T（im）*（（V（im+1）-V（im））/dx+V（im）*（log（A（m+1））-log（A（m）））/dx）;
        pre_rho（m）=rho（im）+d_rho（1m）*dt;        
        pre_V（m）=V（im）+d_V（1m）*dt;
        pre_T（m）=T（im）+d_T（1m）*dt;
    end
    
    % correction step
    d_rho_av=zeros（1nx-1）;
    d_V_av=zeros（1nx-1）;
    d_T_av=zeros（1nx-1）;
    for m=2:nx-1