# matlab 与 python 在科学计算中的区别比较

CFD 专栏收录该内容
1 篇文章 0 订阅

## matlab代码如下(搬运)：

%从亚音速到超音速的转变，非守恒形式
%实现从t时刻到t+dt时刻的转换，rou,V,T
%rou rou2 分别表示rou(i),rou(i+1)
%A,gamma,dt,dx
clear,clc
NUM=31;
C=0.5;%计算时间dt
dx=3/(NUM-1);
step=1000;%时间模拟的步数
Time=0;
gamma=1.4;
x=0:dx:3;
A=1+2.2*(x-1.5).^2;

rou_ar=zeros(step+1,NUM);V_ar=zeros(step+1,NUM);T_ar=zeros(step+1,NUM);
Ma_ar=zeros(step+1,NUM);flux_m=zeros(step+1,NUM);P_ar=zeros(step+1,NUM);
% rou=zeros(1,NUM);V=zeros(1,NUM);T=zeros(1,NUM);
drou_t=zeros(1,NUM);dV_t=zeros(1,NUM);dT_t=zeros(1,NUM);
drou_t2=zeros(1,NUM);dV_t2=zeros(1,NUM);dT_t2=zeros(1,NUM);
drou_t3=zeros(1,NUM);dV_t3=zeros(1,NUM);dT_t3=zeros(1,NUM);
rou_e=zeros(1,NUM);V_e=zeros(1,NUM);T_e=zeros(1,NUM);
rou_r=zeros(1,NUM);V_r=zeros(1,NUM);T_r=zeros(1,NUM);
dt_ar=zeros(1,step);
%初始化
rou=1-0.3146*x;
T=1-0.2314*x;
V=(0.1+1.09*x).*T.^0.5;

%
rou_ar(1,:)=rou;
V_ar(1,:)=V;
T_ar(1,:)=T;
Ma_ar(1,:)=V./sqrt(T);
flux_m(1,:)=rou.*V.*A;
P_ar(1,:)=rou.*T;

for k=1:1:step   %模拟的时间步数

rou=rou_ar(k,:);
V=V_ar(k,:);
T=T_ar(k,:);

%确定本次模拟的dt时间间隔
dt_list=C*dx./(T(2:NUM-1).^0.5+V(2:NUM-1));
dt=min(dt_list);
%
%预测步
for i=2:1:(NUM-1)

drou_t(i)=-rou(i)*(V(i+1)-V(i))/dx-rou(i)*V(i)*(log(A(i+1))-log(A(i)))/dx-V(i)*(rou(i+1)-rou(i))/dx;
dV_t(i)=-V(i)*(V(i+1)-V(i))/dx-1/gamma*((T(i+1)-T(i))/dx+T(i)/rou(i)*((rou(i+1)-rou(i))/dx));
dT_t(i)=-V(i)*(T(i+1)-T(i))/dx-(gamma-1)*T(i)*((V(i+1)-V(i))/dx+V(i)*(log(A(i+1))-log(A(i)))/dx);

rou_e(i)=rou(i)+drou_t(i)*dt;
V_e(i)=V(i)+dV_t(i)*dt;
T_e(i)=T(i)+dT_t(i)*dt;

end
%BC
V_e(1)=2*V_e(2)-V_e(3);
rou_e(1)=1;
T_e(1)=1;

rou_e(NUM)=2*rou_e(NUM-1)-rou_e(NUM-2);
V_e(NUM)=2*V_e(NUM-1)-V_e(NUM-2);
T_e(NUM)=2*T_e(NUM-1)-T_e(NUM-2);
%%%%%%%%%%%%
%预测步结束

%修正步
for i=2:1:(NUM-1)

drou_t2(i)=-rou_e(i)*(V_e(i)-V_e(i-1))/dx-rou_e(i)*V_e(i)*(log(A(i))-log(A(i-1)))/dx-V_e(i)*(rou_e(i)-rou_e(i-1))/dx;
dV_t2(i)=-V_e(i)*(V_e(i)-V_e(i-1))/dx-1/gamma*((T_e(i)-T_e(i-1))/dx+T_e(i)/rou_e(i)*((rou_e(i)-rou_e(i-1))/dx));
dT_t2(i)=-V_e(i)*(T_e(i)-T_e(i-1))/dx-(gamma-1)*T_e(i)*((V_e(i)-V_e(i-1))/dx+V_e(i)*(log(A(i))-log(A(i-1)))/dx);

%对导数求平均
drou_t3(i)=0.5*(drou_t(i)+drou_t2(i));
dV_t3(i)=0.5*(dV_t(i)+dV_t2(i));
dT_t3(i)=0.5*(dT_t(i)+dT_t2(i));

rou_r(i)=rou(i)+drou_t3(i)*dt;
V_r(i)=V(i)+dV_t3(i)*dt;
T_r(i)=T(i)+dT_t3(i)*dt;
end

%BC
V_r(1)=2*V_r(2)-V_r(3);
rou_r(1)=1;
T_r(1)=1;

rou_r(NUM)=2*rou_r(NUM-1)-rou_r(NUM-2);
V_r(NUM)=2*V_r(NUM-1)-V_r(NUM-2);
T_r(NUM)=2*T_r(NUM-1)-T_r(NUM-2);

rou_ar(k+1,:)=rou_r;
V_ar(k+1,:)=V_r;
T_ar(k+1,:)=T_r;
Ma_ar(k+1,:)=V_r./sqrt(T_r);
flux_m(k+1,:)=rou_r.*V_r.*A;
P_ar(k+1,:)=rou_r.*T_r;

dt_ar(k)=dt;
Time=Time+dt;
end  %模拟的时间步数

%理论解的计算
x_theory=0:0.01:3;
A_theory=1+2.2*(x_theory-1.5).^2;
num_x_th=size(x_theory,2);
Ma_theory=zeros(1,num_x_th);
%求解理论上的Ma解
for i=1:1:num_x_th
F=@(ma)(1/ma^2)*((2/(gamma+1))*(1+(gamma-1)/2*ma^2))^((gamma+1)/(gamma-1))-A_theory(i)^2;
if i<=(num_x_th-1)/2+1
Ma_theory(i)=fzero(F,[1e-8,1]);
else
Ma_theory(i)=fzero(F,[1,10]);
end
end

P_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-gamma/(gamma-1));
rou_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1/(gamma-1));
T_theory=(1+((gamma-1)/2)*(Ma_theory(end,:).^2)).^(-1);

%模拟解
P_simu=P_ar(end,:);
rou_simu=rou_ar(end,:);
T_simu=T_ar(end,:);
Ma_simu=Ma_ar(end,:);

%后处理
%针对i节点，显示其在整个模拟过程中的变化。
xn=0:1:step;
i=(NUM-1)/2+1;%设置所要观察的节点
rou_i=rou_ar(:,i);
V_i=V_ar(:,i);
T_i=T_ar(:,i);
Ma_i=Ma_ar(:,i);

subplot(2,2,1)
plot(xn,rou_i);
title('密度');

subplot(2,2,2)
plot(xn,V_i);
title('V');

subplot(2,2,3)
plot(xn,T_i);
title('T');

subplot(2,2,4)
plot(xn,Ma_i);
title('Ma');

%质量流量的变化
figure(2)
flux_list=[0,50,150,200,700]+1;%注：所研究的最大700，不能超过step
plot(x,flux_m(flux_list(1),:),x,flux_m(flux_list(2),:),x,flux_m(flux_list(3),:),x,flux_m(flux_list(4),:),x,flux_m(flux_list(5),:));
legend('0dt','50dt','150dt','200dt','700dt');
title('质量流量');
%稳态的理论解和模拟解的对比
%理论解

%画图对比
figure(3)
subplot(2,2,1)
plot(x_theory,P_theory,'-',x,P_simu,'+');
legend('理论解','模拟解');
title('Pressure comparision');

subplot(2,2,2)
plot(x_theory,rou_theory,'-',x,rou_simu,'+');
legend('理论解','模拟解');
title('Density comparision');

subplot(2,2,3)
plot(x_theory,T_theory,'-',x,T_simu,'+');
legend('理论解','模拟解');
title('Temperature comparision');

subplot(2,2,4)
plot(x_theory,Ma_theory,'-',x,Ma_simu,'+');
legend('理论解','模拟解');
title('Mach comparision');


## 与之对应的python3 源码如下（原创）：

import numpy
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
for i in range(10):
print(i)

NUM = 31
C = 0.5
dx = 3/(NUM-1)
step = 1000
Time = 0
gamma = 1.4
x = numpy.arange(0, 3.1, dx)
A = 1 + 2.2 * pow((x-1.5),2)

rou_ar = numpy.zeros([step+1, NUM])
V_ar = numpy.zeros([step+1, NUM])
T_ar = numpy.zeros([step+1, NUM])
Ma_ar = numpy.zeros([step+1,NUM])
flux_m =numpy.zeros([step+1,NUM])
P_ar = numpy.zeros([step+1,NUM])

drou_t = numpy.zeros([1,NUM])
dV_t =numpy.zeros([1,NUM])
dT_t = numpy.zeros([1,NUM])
drou_t2 = numpy.zeros([1,NUM])
dV_t2 = numpy.zeros([1,NUM])
dT_t2 = numpy.zeros([1,NUM])
drou_t3 = numpy.zeros([1,NUM])
dV_t3 = numpy.zeros([1,NUM])
dT_t3 = numpy.zeros([1,NUM])
rou_e =numpy.zeros([1,NUM])
V_e = numpy.zeros([1,NUM])
T_e = numpy.zeros([1,NUM])
rou_r = numpy.zeros([1,NUM])
V_r = numpy.zeros([1,NUM])
T_r = numpy.zeros([1,NUM])
dt_ar=numpy.zeros([1,step])
print(drou_t[0,1])

rou = 1 - 0.3146 * x
T = 1 - 0.2314 * x
V = (0.1 + 1.09 * x) * numpy.sqrt(T)
print(V[0])
for i in range(NUM):
rou_ar[0,i] = rou[i]
V_ar[0,i]=V[i]
T_ar[0,i]=T[i]
Ma_ar[0,i]=V[i]/numpy.sqrt(T[i])
flux_m[0,i]=rou[i]*V[i]*A[i]
P_ar[0,i]=rou[i]*T[i]

for k in range(step):
for j in range(NUM):
rou[j] = rou_ar[k,j]
V[j] = V_ar[k,j]
T[j] = T_ar[k,j]
print(T)
print(T[0])
print(T[1:5])
t = numpy.sqrt(T[1:])
v = numpy.sqrt(V[1:])
dt_list = C * dx/(numpy.sqrt(T[1:])+ V[1:])
# print(dt_list)
dt = min(dt_list)

for i in range(1,NUM-1,1):
drou_t[0,i] = -rou[i] * (V[i+1] - V[i]) / dx - rou[i] * V[i] * (numpy.log(A[ i + 1]) - numpy.log(A[i])) / dx - V[i] * (rou[i + 1] - rou[i]) / dx
dV_t[0,i] = -V[i] * (V[i+1] - V[i]) / dx - 1/gamma * (((T[i+1]-T[i])/dx) + T[i]/rou[i] * ((rou[i+1] - rou[i]) / dx))
dT_t[0,i] = -V[i] * (T[i + 1] - T[i]) / dx - (gamma - 1) * T[i] * ((V[i + 1] - V[i]) / dx + V[i] * (numpy.log(A[i + 1]) - numpy.log(A[i])) / dx)
rou_e[0,i] = rou[i] + drou_t[0,i] * dt
V_e[0,i] = V[i] + dV_t[0,i] * dt
T_e[0,i] = T[i] + dT_t[0,i] * dt

V_e[0,0] = 2 * V_e[0,1] - V_e[0,2]
rou_e[0,0] = 1
T_e[0,0] = 1
rou_e[0,NUM-1] = 2 * rou_e[0,NUM - 2] - rou_e[0,NUM - 3]
V_e[0,NUM-1] = 2 * V_e[0,NUM - 2] - V_e[0,NUM - 3]
T_e[0,NUM-1] = 2 * T_e[0,NUM - 2] - T_e[0,NUM - 3]

for i in range(1, NUM-1, 1):
drou_t2[0,i] = -rou_e[0,i] * (V_e[0,i] - V_e[0,i - 1]) / dx - rou_e[0,i] * V_e[0,i] * (numpy.log(A[i]) - numpy.log(A[i - 1])) / dx - V_e[0,i] * (rou_e[0,i] - rou_e[0,i - 1]) / dx
dV_t2[0,i] = -V_e[0,i] * (V_e[0,i] - V_e[0,i - 1]) / dx - 1 / gamma * ((T_e[0,i] - T_e[0,i - 1]) / dx + T_e[0,i] / rou_e[0,i] * ((rou_e[0,i] - rou_e[0,i - 1]) / dx))
dT_t2[0,i] = -V_e[0,i] * (T_e[0,i] - T_e[0,i - 1]) / dx - (gamma - 1) * T_e[0,i] * ((V_e[0,i] - V_e[0,i - 1]) / dx + V_e[0,i] * (numpy.log(A[i]) - numpy.log(A[i - 1])) / dx)

drou_t3[0,i] = 0.5 * (drou_t[0,i] + drou_t2[0,i])
dV_t3[0,i] = 0.5 * (dV_t[0,i] + dV_t2[0,i])
dT_t3[0,i] = 0.5 * (dT_t[0,i] + dT_t2[0,i])

rou_r[0,i] = rou[i] + drou_t3[0,i] * dt
V_r[0,i] = V[i] + dV_t3[0,i] * dt
T_r[0,i] = T[i] + dT_t3[0,i] * dt

V_r[0,0] = 2 * V_r[0,1] - V_r[0,2]
rou_r[0,0] = 1
T_r[0,0] = 1
rou_r[0,NUM - 1] = 2 * rou_r[0,NUM - 2] - rou_r[0,NUM - 3]
V_r[0,NUM - 1] = 2 * V_r[0,NUM - 2] - V_r[0,NUM - 3]
T_r[0,NUM - 1] = 2 * T_r[0,NUM - 2] - T_r[0,NUM - 3]

for m in range(NUM):
rou_ar[k+1, m] = rou_r[0,m]
V_ar[k + 1,m]=V_r[0,m]
T_ar[k + 1,m]=T_r[0,m]
Ma_ar[k + 1,m]=V_r[0,m]/ numpy.sqrt(T_r[0,m])
flux_m[k + 1,m]=rou_r[0,m]* V_r[0,m]* A[m]
P_ar[k + 1,m]=rou_r[0,m]* T_r[0,m]
dt_ar[0,k] = dt
Time = Time + dt

x_theory = numpy.arange(0,3.01,0.01)
A_theory = 1 + 2.2 * pow((x_theory - 1.5),2)
# num_x_th = x_theory.shape()

x_theory = numpy.arange(0,3.01,0.01)
A_theory = 1 + 2.2 * pow((x_theory - 1.5),2)
num_x_th = numpy.shape(x_theory)
print(num_x_th[0])
a = num_x_th[0]
Ma_theory = numpy.zeros([1,num_x_th[0]])
gamma = 1.4
print(A_theory)
# fsolve(f,[1,2])
ma = 0.0978
print("ma",Ma_theory[0,1])
i=0
f=(1 / ma ** 2) * pow((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * ma ** 2), ((gamma + 1) / (gamma - 1)))- A_theory[i]**2
print("f=",f)
print(A_theory[i]**2)
print(pow(A_theory[i],2))
def f(ma,i):
return (1 / ma ** 2) * pow(((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * ma ** 2)), ((gamma + 1) / (gamma - 1)))- A_theory[i]**2
#
# print(fsolve(f,[-2,2]))
#
#
bb = numpy.array([1e-8,1])
for i in range(301):
if i <=(num_x_th[0] - 1)/2 + 1:
Ma_theory[0,i]=fsolve(f,0.01,i)
aa=Ma_theory[0,i]
print(i,aa)
else:
Ma_theory[0,i] = fsolve(f, 1.5, i)
print(i, aa)

P_theory=pow((1+((gamma-1)/2)*(Ma_theory[0,:]**2)),(-gamma/(gamma-1)))
rou_theory=pow((1+((gamma-1)/2)*(Ma_theory[0:]**2)),(-1/(gamma-1)))
T_theory=(1+((gamma-1)/2)/(Ma_theory[0,:]**2))
print("P", P_theory)

P_simu=P_ar[-1,:]
rou_simu=rou_ar[-1,:]
T_simu=T_ar[-1,:]
Ma_simu=Ma_ar[-1,:]
print(P_simu-P_ar[300,:])

xn = numpy.arange(0,step+1,1)
i = (NUM - 1) / 2 + 1
ii= int(i)
print(rou_ar[:,ii])

rou_i = rou_ar[:,ii]
print(rou_ar[:,ii])
V_i = V_ar[:, ii]
T_i = T_ar[:, ii]
Ma_i = Ma_ar[:,ii]

plt.subplot(2,2,1)
plt.plot(xn.T, rou_i)
plt.title("密度")

plt.subplot(2,2,2)
plt.plot(xn.T, V_i)
plt.title("V")

plt.subplot(2,2,3)
plt.plot(xn.T, T_i)
plt.title("T")

plt.subplot(2,2,4)
plt.plot(xn.T, Ma_i)
plt.title("Ma")

plt.show()


ps:python源码中省略了后处理的过程。

python的缺点在于调试、运行速度比较慢，对于我这样的小白，不知道如何利用高级python功能的话极易造成程序冗长，程序运行相对缓慢。

python最致命的缺点在于有些错误位置找不出来，因而修改不便。

python要注意：
x=range(10),那么x是一维数组，按索引查看元素过程中，一定是x[0], x[1], x[2];
y=numpy.zeros(1,2), 生成的y是1行2列的二维数组，引用一定是 y[0,0], y[0,1]
x=30, y=10, a=x/y 在python中可能得出的a 是浮点型，用于数组的索引时程序报错但不会和你说具体原因的，千万要转成int

matlab中的fsolve函数可以求指定区间上函数的零点，fsolve(myfun,[0,5])–求myfun函数在区间[0,5]范围内的零点。

fsolve(myfun,start_point),计算从一个点开始最接近改点的函数零点。如果在区间内存在多个零点，那么python中的fsolve如何使用，还请大神们指教。

• 2
点赞
• 2
评论
• 2
收藏
• 一键三连
• 扫一扫，分享海报

04-27 727
04-10 2122
12-16 3436
08-03