使用python改写matlab程序

准备将这段matlab程序用pycharm来写

Riemann问题精确解及程序实现_梅冠华的博客-CSDN博客_一维riemann问题的求解一维无粘流动初始间断的演化问题，激波管问题（Sod激波管问题），由于该问题属于典型的间断问题，且有精确解存在，故广泛用于对比验证CFD中离散格式、数值方法的准确性，意义嘛还是蛮重要的。1 问题描述 2 求解方法 3 解法流程 4 程序实现 5 计算结果

https://blog.csdn.net/meiguanhua/article/details/104454790

我对python中数组、列表的使用并不明确，总是出现各种报错，不知道如何修改可以达到matlab中的效果。

  File "D:\Pycharm\learnpython\Riemann1.0.py", line 232, in 
    xSpn = np.sum(xSpnLftWavL, xSpnLftExp, xSpnMidL, xSpnMidR, xSpnRgtExp, xSpnRgtWavR)
TypeError: Cannot construct a dtype from an array

问题相关代码，请勿粘贴截

import sympy as sm
import numpy as np
import math
import matplotlib.pyplot as plt

tNow = 0.14
gamma = 1.4
U1 = 0
R1 = 1
P1 = 1
U2 = 0
R2 = 0.125
P2 = 0.1
C1 = sm.sqrt(gamma*P1/R1)
C2 = sm.sqrt(gamma*P2/R2)

# 定义函数F(p*) = f(p*, p1, rho1) + f(p*, p2, rho2)
def function1(Ps, Pi, Ri):
    gamma = 1.4
    Ci = sm.sqrt(gamma*Pi/Ri)
    if(Ps>Pi):
        temp = ((gamma+1)/(2*gamma))*(Ps/Pi)+((gamma-1)/(2*gamma))
        f = (Ps-Pi)/((Ri*Ci)*sm.sqrt(temp))
    elif(Pstemp = (Ps/Pi)**((gamma-1)/(2*gamma))-1
        f = (2*Ci)/(gamma-1)*temp
    else:
        f = 0
    return f
# 判断情况的函数
def function2(Fp0, Fp1, Fp2, U1, U2, P1, P2):
    U1_U2 = U1-U2
    idxCase = 0
    if (P2>=P1):
        if ((U1_U2)>=Fp2):
            idxCase = 1
        elif(Fp1<=(U1_U2)idxCase = 3
        elif(Fp0<=(U1_U2)idxCase = 4
        elif(U1_U2idxCase = 5
    elif (P2if ((U1_U2)>=Fp1):
            idxCase = 1
        elif(Fp2<=(U1_U2)idxCase = 2
        elif(Fp0<=(U1_U2)idxCase = 4
        elif(U1_U2idxCase = 5
    return idxCase
# 二分法求解中间区域压力和速度
def function3(U1, P1, R1, U2, P2, R2):
    global Ps
    PSpn = range(0, 10000)  # 压力的范围
    U1_U2 = U1-U2
    for iter in range(0, 1000):
        Ps = np.mean(PSpn)
        Fps = function1(Ps, P1, R1)+function1(Ps, P2, R2)
        if (U1_U2>Fps):
            PSpn = np.arange(Ps, 10000)
        elif (U1_U2PSpn = np.arange(0, Ps)
#        if ((PSpn[1]-PSpn[0])<(1e-10)):
#            Ps = np.mean(PSpn)
    Us = 0.5*(U1+U2+function1(Ps, P2, R2)-function1(Ps, P1, R1))
    return Ps, Us
# 左激波情况下计算激波运行速度Z1和中间区域左侧密度R1s
def function4(P1, U1, R1, Ps, Us, gamma):
    Z1 = U1-(Ps-P1)/(R1*(U1-Us))
    R1s = R1*(U1-Z1)/(Us-Z1)
    return Z1, R1s
# 左膨胀波情况下计算波头波尾运行速度Z1Head, Z1Tail
def function5(P1, U1, R1, Ps, Us, gamma, idxCase):
    C1 = sm.sqrt(gamma * P1 / R1)
    Z1Head = U1 - C1
    if (idxCase!=5):
        R1s = R1 * (Ps / P1) ** (1 / gamma)
        C1s = sm.sqrt(gamma * Ps / R1s)
        Z1Tail = Us - C1s
    else:
        R1s = 0
        Z1Tail = U1-2*C1/(gamma-1)
    return R1s, Z1Head, Z1Tail
# 右激波情况下计算激波运行速度Z2和中间区域右侧密度R2s
def function6(P2, U2, R2, Ps, Us):
    Z2 = U2-(Ps-P2)/(R2*(U2-Us))
    R2s = R2*(U2-Z2)/(Us-Z2)
    return Z2, R2s
# 右膨胀波情况下计算波头波尾运行速度Z2Head, Z2Tail
def function7(P2, U2, R2, Ps, Us, gamma, idxCase):
    C2 = sm.sqrt(gamma*P2/R2)
    Z2Head = U2+C2
    if (idxCase!=5):
        R2s = R2*(Ps/P2)**(1/gamma)
        C2s = sm.sqrt(gamma*Ps/R2s)
        Z2Tail = Us+C2s
    else:
        R2s = 0
        Z2Tail = U2+2*C2/(gamma-1)
    return R2s, Z2Head, Z2Tail
# 左膨胀波内部状态值u , p , ρ
def function8(Z1Head, Z1Tail, tNow, gamma, U1, C1, P1):
    xDwn = Z1Head*tNow
    xUp = Z1Tail*tNow
    xDlt = (xUp-xDwn)/100
    xSpnLftExp = np.arange(xDwn, xUp, xDlt)  # np.arange(起点, 终点, 步长)
    cSpnLftExp = (gamma-1)/(gamma+1)*(U1-xSpnLftExp/tNow)+2/(gamma+1)*C1
    uSpnLftExp = xSpnLftExp/tNow+cSpnLftExp
    pSpnLftExp = P1*(cSpnLftExp/C1)**(2*gamma/(gamma-1))
    rSpnLftExp = gamma*pSpnLftExp/cSpnLftExp**2
    return xSpnLftExp, uSpnLftExp, pSpnLftExp, rSpnLftExp
# 右膨胀波内部状态值u , p , ρ
def function9(Z2Tail, Z2Head, tNow, gamma, U2, C2, P2):
    xDwn = Z2Tail*tNow
    xUp = Z2Head*tNow
    xDlt = (xUp-xDwn)/100
    xSpnRgtExp = np.arange(xDwn, xUp, xDlt)
    cSpnRgtExp = -(gamma-1)/(gamma+1)*(U2-xSpnRgtExp/tNow)+2/(gamma+1)*C2
    uSpnRgtExp = xSpnRgtExp/tNow-cSpnRgtExp
    pSpnRgtExp = P2*(cSpnRgtExp/C2)**(2*gamma/(gamma-1))
    rSpnRgtExp = gamma*pSpnRgtExp/cSpnRgtExp**2
    return xSpnRgtExp, uSpnRgtExp, pSpnRgtExp, rSpnRgtExp
# 计算左侧区域范围的函数
def function10(Z1, P1, U1, R1, tNow):
    lftShockX = Z1*tNow
    xSpnLftWavL = np.arange(min(-1, int(lftShockX*2)), lftShockX)
    pSpnLftWavL = np.arange(P1, P1)
    uSpnLftWavL = np.arange(U1, U1)
    rSpnLftWavL = np.arange(R1, R1)
    return xSpnLftWavL, pSpnLftWavL, uSpnLftWavL, rSpnLftWavL
# 计算中间区域位于间断左侧范围的函数
def function11(Z1, Ps, Us, RsL, tNow):
    lftShockX = Z1*tNow
    midDisconX = Us*tNow
    xSpnMidL = np.arange(lftShockX, midDisconX)
    pSpnMidL = np.arange(Ps, Ps)
    uSpnMidL = np.arange(Us, Us)
    rSpnMidL = np.arange(RsL, RsL)
    return xSpnMidL, pSpnMidL, uSpnMidL, rSpnMidL
# 计算中间区域位于间断右侧范围的函数
def function12(Z2, Ps, Us, RsR, tNow):
    rgtShockX = Z2*tNow
    midDisconX = Us*tNow
    xSpnMidR = np.arange(midDisconX, rgtShockX)
    pSpnMidR = np.arange(Ps, Ps)
    uSpnMidR = np.arange(Us, Us)
    rSpnMidR = np.arange(RsR, RsR)
    return xSpnMidR, pSpnMidR, uSpnMidR, rSpnMidR
# 计算右侧区域范围的函数
def function13(Z2, P2, U2, R2, tNow):
    rgtShockX = Z2*tNow
    xSpnRgtWavR = np.arange(rgtShockX, max(1, math.ceil(rgtShockX*2)))
    pSpnRgtWavR = np.arange(P2, P2)
    uSpnRgtWavR = np.arange(U2, U2)
    rSpnRgtWavR = np.arange(R2, R2)
    return xSpnRgtWavR, pSpnRgtWavR, uSpnRgtWavR, rSpnRgtWavR

Fp0 = function1(0, P1, R1)+function1(0, P2, R2)
Fp1 = function1(P1, P1, R1)+function1(P1, P2, R2)
Fp2 = function1(P2, P1, R1)+function1(P2, P2, R2)

idxCase = function2(Fp0, Fp1, Fp2, U1, U2, P1, P2)

Ps, Us = function3(U1, P1, R1, U2, P2, R2)

global Z1, Z2, Z1Head, Z1Tail, Z2Tail, Z2Head

# 中心区接触间断两侧的密度,以及左右波的传播速度
# 左波为激波的情况（情况1和3）
if ((idxCase == 1)or(idxCase == 3)):
    Z1, RLs = function4(P1, U1, R1, Ps, Us, gamma)
# 左波为膨胀波（稀疏波）的情况（情况2, 4, 5）
else:
    RLs, Z1Head, Z1Tail = function5(P1, U1, R1, Ps, Us, gamma, idxCase)

# 右波为激波的情况（情况1和2）
if ((idxCase == 1)or(idxCase == 2)):
    Z2, R2s = function6(P2, U2, R2, Ps, Us)
# 右波为膨胀波（稀疏波）的情况（情况3, 4, 5）
else:
    R2s, Z2Head, Z2Tail = function7(P2, U2, R2, Ps, Us, gamma, idxCase)

# 计算稀疏波区域的值
# 左稀疏波（情况2, 4, 5）
if ((idxCase == 2)or(idxCase == 4)):
    xSpnLftExp, uSpnLftExp, pSpnLftExp, rSpnLftExp = function8(Z1Head, Z1Tail, tNow, gamma, U1, C1, P1)
else:
    xSpnLftExp = []
    uSpnLftExp = []
    pSpnLftExp = []
    rSpnLftExp = []

# 右稀疏波（情况3, 4, 5）
if ((idxCase == 3)or(idxCase == 4)):
    xSpnRgtExp, uSpnRgtExp, pSpnRgtExp, rSpnRgtExp = (function9(Z2Tail, Z2Head, tNow, gamma, U2, C2, P2))[0]
else:
    xSpnRgtExp = []
    uSpnRgtExp = []
    pSpnRgtExp = []
    rSpnRgtExp = []

# 左激波或左膨胀波左侧区域
# 左波为激波的情况（情况1和3）
if ((idxCase == 1)or(idxCase == 3)):
    xSpnLftWavL, pSpnLftWavL, uSpnLftWavL, rSpnLftWavL = function10(Z1, P1, U1, R1, tNow)
else:
    xSpnLftWavL, pSpnLftWavL, uSpnLftWavL, rSpnLftWavL = function10(Z1Head, P1, U1, R1, tNow)

# 中间区域（接触间断左侧区域）
# 左波为激波的情况（情况1和3）
if ((idxCase == 1)or(idxCase == 3)):
    xSpnMidL, pSpnMidL, uSpnMidL, rSpnMidL = function11(Z1, Ps, Us, RLs, tNow)
else:
    xSpnMidL, pSpnMidL, uSpnMidL, rSpnMidL = function11(Z1Tail, Ps, Us, RLs, tNow)

# 中间区域（接触间断面右侧）
# 右波为激波的情况（情况1和2)
if ((idxCase == 1)or(idxCase == 2)):
    xSpnMidR, pSpnMidR, uSpnMidR, rSpnMidR = function12(Z2, Ps, Us, R2s, tNow)
else:
    xSpnMidR, pSpnMidR, uSpnMidR, rSpnMidR = function12(Z2Tail, Ps, Us, R2s, tNow)

# 右激波或膨胀波右侧区域
# 右波为激波的情况（情况1和2）
if ((idxCase == 1)or(idxCase == 2)):
    xSpnRgtWavR, pSpnRgtWavR, uSpnRgtWavR, rSpnRgtWavR = function13(Z2, P2, U2, R2, tNow)
else:
    xSpnRgtWavR, pSpnRgtWavR, uSpnRgtWavR, rSpnRgtWavR = function13(Z2Head, P2, U2, R2, tNow)

xSpn = np.sum(xSpnLftWavL, xSpnLftExp, xSpnMidL, xSpnMidR, xSpnRgtExp, xSpnRgtWavR)
uSpn = np.sum(uSpnLftWavL, uSpnLftExp, uSpnMidL, uSpnMidR, uSpnRgtExp, uSpnRgtWavR)
pSpn = np.sum(pSpnLftWavL, pSpnLftExp, pSpnMidL, pSpnMidR, pSpnRgtExp, pSpnRgtWavR)
rSpn = np.sum(rSpnLftWavL, rSpnLftExp, rSpnMidL, rSpnMidR, rSpnRgtExp, rSpnRgtWavR)

# 绘制图形
plt.figure(99)
plt.plot(xSpn, pSpn, linewidth=2)
plt.plot(xSpn, uSpn, linewidth=2)
plt.plot(xSpn, rSpn, linewidth=2)
plt.legend('p', 'u', 'rho')
plt.xlabel('x')
plt.ylabel('p, u, rho')
plt.title('Riemann Problem Accurate Solution')
plt.show()

把他改成中括号试试？