使用 odeint 求解 odes

Question

我最近尝试用 python 练习颂歌。下面的代码是关于求解三个相互关联的方程式。也就是说，第一个值需要是下一个等式的值，依此类推。

这是我的代码，但是 odeint 有问题并打印警告 'setting an array element with a sequence.'

迷茫了好久，谁能看出错误吗？

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

rho=1000.0  
A=10.0
g=9.8   
Gf=1.0 
fi=10.0
f0=2.0
n=101
f_=0.0
h_=1.63265306e-05 
t=np.linspace(0,60,n)
cmax=2.0

def model(z,t,v): 

    level1=z[0]
    level2=z[1]
    level3=z[2]
    v11=v[0]
    v22=v[1]
    v33=v[2]

    dlevel1_dt=(fi-f0-((cmax*v11)*((rho*g*level1/Gf)**(1/2))))/A
    f1=f_+(cmax*v1*((rho*g/Gf/h_)**2))/2*(level1-h_)

    dlevel2_dt=(f1-((cmax*v22)*((rho*g*level2/Gf)**(1/2))))/A
    f2=f_+(cmax*v2*((rho*g/Gf/h_)**2))/2*(level2-h_) 

    dlevel3_dt=(f2-((cmax*v33)*((rho*g*level3/Gf)**(1/2))))/A
    f3=f_+(cmax*v3*((rho*g/Gf/h_)**2))/2*(level3-h_) 

    dlevel_dt=[dlevel1_dt,dlevel2_dt,dlevel3_dt] 
    return dlevel_dt


level1=np.empty_like(t)
level2=np.empty_like(t)
level3=np.empty_like(t)

level0=[0,0,0]
level1[0]=level0[0]
level2[0]=level0[1]
level3[0]=level0[2]

v1=np.zeros(n)
v1[5:20]=10 
v2=np.zeros(n)
v2[30:50]=5 
v3=np.zeros(n)
v3[65:85]=2 
v=[v1,v2,v3]

for i in range(1,n): 
    tspan=[t[i-1],t[i]]    
    z=odeint(model,level0,tspan,args=(v[i],))
    level1[i]=z[1][0]
    level2[i]=z[1][1]
    level3[i]=z[1][2]   
    level0=z[1]

plt.figure()         
plt.plot(t,level1,t,level2,t,level3)
plt.show()

结果：

z=odeint(model,level0,tspan,args=(v[i],)) line 244, in odeint
int(bool(tfirst)) ValueError: setting an array element with a sequence.

Answer 1

完整的错误信息

Traceback (most recent call last):
  File "<ipython-input-1-49ffe15670c0>", line 56, in <module>
    z=odeint(model,level0,tspan,args=(v[i],))
  File "/usr/local/lib/python3.8/dist-packages/scipy/integrate/odepack.py", line 241, in odeint
    output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
ValueError: setting an array element with a sequence. 
The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

通过对 model 的示例调用：

In [4]: res = model(level0,0,v[0])
In [5]: len(res)
Out[5]: 3
In [6]: res[0]
Out[6]: 0.8
In [7]: res[1].shape
Out[7]: (101,)
In [8]: res[2].shape
Out[8]: (101,)

return 是一个 3 元素列表。第一个元素是标量，但其他 2 个是数组。 odeint 无法使用那种 model。

对于dlevel1_dt，fi是一个数字。

对于dlevel2_dt，f1是从v1推导出来的，是一个(101,)形数组。因此 dlevel2_dt 是一个数组。

odeint 试图将 model 的结果转换为 float dtype 数组：

In [13]: np.array(res, float)
Traceback (most recent call last):
  File "<ipython-input-13-703f82987b55>", line 1, in <module>
    np.array(res, float)
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

总之，model 应该 return 一个在形状上匹配 z 的浮点 numpy 数组。

使用 odeint 求解 odes

using odeint to solve odes

python

numpy

scipy

odeint

完整的错误信息