拟一维喷管流动的数值解MATLAB程序

clear
clc
tic
%-----Code By WangNianhua 2013.10.17
% 亚音速-超音速等熵喷管，非守恒形式方程MacCormack方法求解，教材P195-222
%---------------第一步，喷管形状参数和初始条件--------
L=3;N=31;dx=L/（N-1）;C=0.5;
x=0:dx:L; 
A=1+2.2*（x-1.5）.^2;
subplot（221）plot（xA‘-‘）;
hold on;
plot（x-A‘-‘）;
plot（xzeros（1N）‘.-‘）
title（‘喷管外形（网格点分布如图）‘）
hold off;
Rou=1-0.3146*x;
T=1-0.2314*x;
V=（0.1+1.09*x）.*sqrt（T）;
for n=1:1400
    Y1=Rou（16）;
%-------------第二步，预估步------------------------
for i=2:30
    P_rou_t（i）=-V（i）*（（Rou（i+1）-Rou（i））/dx）-Rou（i）*（（V（i+1）-V（i））/dx）-Rou（i）*V（i）*（（log（A（i+1））-log（A（i）））/dx）;
    P_v_t（i）=-V（i）*（（V（i+1）-V（i））/dx）-（（T（i+1）-T（i））/dx+（（Rou（i+1）-Rou（i））/dx）*T（i）/Rou（i））*1/1.4;
    P_T_t（i）=-V（i）*（（T（i+1）-T（i））/dx）-0.4*T（i）*（（（V（i+1）-V（i））/dx）+V（i）*（（log（A（i+1））-log（A（i）））/dx））;
end
%-------------第三步，求Δt求rou V T的预估值-----
dt=C*（dx./（V（2:30）+sqrt（T（2:30））））;
dt=min（dt）;
Rou1（2:30）=Rou（2:30）+P_rou_t（2:30）.*dt;
V1（2:30）=V（2:30）+P_v_t（2:30）.*dt;
T1（2:30）=T（2:30）+P_T_t（2:30）.*dt;
V1（1）=V（1）;T1（1）=T（1）;Rou1（1）=Rou（1）;
%--------------第四步，校正步--------------
for i=2:30
    P_rou_t_2（i）=-V1（i）*（（Rou1（i）-Rou1（i-1））/dx）-Rou1（i）*（（V1（i）-V1（i-1））/dx）-Rou1（i）*V1（i）