我是微积分和 Python/Scipy 的新手，所以如果这个问题太愚蠢，我深表歉意。我正在尝试模拟两个压力容器之间的流动。假设我们有两个点和它们之间的链接，就像这样

[$Vc_1$,$P_1$]----($A$)----[$Vc_2$,$P_2$]

$Vc_1$,$Vc_2$是节点（容器）的常数体积和$P_1$,$P_2$分别是节点的不同压力。

我最终在下面写了不同的问题。没关系物理意义，我只想让数学正确。

$\frac{\mathrm{dP_1} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_1 \cdot B}$

$\frac{\mathrm{dP_2} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_2 \cdot B}$

这里$B$是可压缩性。

$\frac{\mathrm{dVa} }{\mathrm{d} t} = A \cdot K \cdot \sqrt {P_1 - P_2}$

$\frac{\mathrm{dVa} }{\mathrm{d} t}$是节点之间的“流量”或附加体积的变化量。$K$是一些常数系数和$A$是链接“吞吐量”。($P_1-P_2$) 可以更改符号，因此我已在软件中对此进行了调整。

下面是我编写的用来评估它的 Python 程序。

#!/usr/bin/env python

导入数学
从 scipy.integrate 导入 o​​deint
从时间进口时间
导入 numpy

B_compressibility = 0.0000033 #水的可压缩性
K = 0.747871759938 #系数

Vc_1 = 20
Vc_2 = 50
A = 0.01
P_1 = 4000
P_2 = 2000

def deriv（状态，t）：    
    _P_1 = 状态[0]
    _P_2 = 状态[2]   
    diff_P = _P_1 - _P_2
    flow_direction = math.copysign(1, diff_P)
    dVa = flow_direction * A * K * math.sqrt(abs(diff_P))
    dP_1 = -(dVa/Vc_1)/B_compressibility
    dP_2 = (dVa/Vc_2)/B_compressibility
    #print 'IN'，状态
    #print 'OUT', [dP_1, -dVa, dP_2, dVa]
    返回 [dP_1，-dVa，dP_2，dVa]

如果 __name__ == '__main__'：
    Va_1 = Vc_1 * P_1 * B_compressibility
    Va_2 = Vc_2 * P_2 * B_compressibility
    odeIterations = 10
    时间段 = numpy.linspace(0.0,1.0, odeIterations)    
    初始状态 = [P_1，Va_1，P_2，Va_2]
    t0 = 时间（）
    state_array = odeint（衍生，初始状态，时间段）
    t1 = 时间（）
    打印“运行时 %fs”%(t1-t0)        
    打印状态数组
    P_1，Va_1，P_2，Va_2 = state_array[odeIterations-1]

下面是程序的输出

在此调用上完成的工作过多（可能是错误的 Dfun 类型）。
以 full_output = 1 运行以获取定量信息。
运行时间 0.041000s
[[ 4.00000000e+003 2.64000000e-001 2.00000000e+003 3.30000000e-001]
 [ 3.49242034e+003 2.30499743e-001 2.20303186e+003 3.63500257e-001]
 [ 3.09580400e+003 2.04323064e-001 2.36167840e+003 3.89676936e-001]
 [ 2.81015098e+003 1.85469965e-001 2.47593961e+003 4.08530035e-001]
 [ 2.63546127e+003 1.73940444e-001 2.54581549e+003 4.20059556e-001]
 [ 2.57173487e+003 1.69734501e-001 2.57130605e+003 4.24265499e-001]
 [ 2.57142857e+003 1.69714286e-001 2.57142857e+003 4.24285714e-001]
 [1.83357137e-299 1.80790662e-299 1.83145695e-299 1.87152166e-299]
 [1.83276935e-299 1.80681296e-299 1.83182150e-299 1.87141230e-299]
 [1.83379011e-299 1.80746916e-299 1.83320682e-299 1.83229543e-299]]
 lsoda-- 在当前 t (=r1), mxstep (=i1) 步   
       在进行宣传之前接听此电话     
      在上述消息中，I1 = 500
      在上述消息中，R1 = 0.6110315150411E+00

odeint给出正确的结果直到第 7 行，然后出现严重错误。我在谷歌搜索过，看起来我不是唯一一个与 scipy 斗争的人。每个人都建议增加 mxstep 但这并不能解决我的问题。此外，它显着减慢了该方法。有人建议降低准确性，但我不知道该怎么做。如果有帮助，降低准确性是可以的，因为我不需要odeint. 点后的几个数字对我来说绰绰有余。另外我只需要最终值，所以实际上我想减少步骤数。任何帮助是极大的赞赏！