【发布时间】：2020-11-01 10:56:07
【问题描述】：
我正在尝试解决以下一对 ODE：
- dq/dt = k
- dk/dt = 1/L (E(t) - R(i) k - 1/C q)
这个想法是求解 q，然后绘制 i(t) = dq/dt。完整代码如下：
import timeit
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
start = timeit.default_timer()
#Data:
t = np.linspace(0, 10, 100)
w1 = 350 #350 rad/s
w2 = 1e3 #1000 rad/s
phi1 = 150
phi2 = 30
E = np.cos(t)
R1 = 10e3 #10 kOhm
R2 = 3.3e3 #3k3 Ohm
L = 10e-3 #10 mHy
C = 1.56e-6 #1.56 uF
#Create model -> x = [q , k]
def model(x , t):
q = x[0]
k = x[1]
dq_dt = x[1]
i = dq_dt
R = R1 * i + R2 * i**3
dk_dt = 1 / L * (E - R * k - 1/C * q)
dx_dt = np.array([dq_dt , dk_dt])
return dx_dt
#Init. cond.:
x0 = np.array([0 , 0])
#Solver ODE:
x = odeint(model, x0, t)
q = x[: , 0]
i = np.diff(q)
vr = (R1 * i + R2 * i**3) * i
vl = L * np.diff(i)
vc = 1/C * q
#Plot:
plt.style.use("bmh")
fig,(ax1,ax2) = plt.subplots(nrows = 2, ncols = 1, sharex = True)
ax1.plot(t, E, "r-", linewidth = 2)
ax2.plot(t, np.append(i,0), "b-", linewidth = 2)
ax1.grid(True)
ax1.set_title("E(t)")
ax1.set_ylabel("E(t) [V]")
ax2.grid(True)
ax2.set_title("i(t)")
ax2.set_ylabel("i(t) [A]")
ax2.set_xlabel("t [s]")
plt.tight_layout()
plt.show()
end = timeit.default_timer()
print ("\n\n\nSIM TIME = ", end - start, " s")
当我像这样运行代码时，出现以下错误：
File "...", line 73, in <module>
x = odeint(model, x0, t)
File ".../python3.7/site-packages/scipy/integrate/odepack.py", line 244, in odeint
int(bool(tfirst)))
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'
问题在于dx_dt = np.array([dq_dt , dk_dt]) 的返回。由于E 是一个包含100 个元素的数组，dk_dt = dx_dt[1] 也将是一个包含100 个元素的数组，而dq_dt = dx_dt[0] 只是一个值。如果我这样做 dx_dt = np.array([dq_dt , dk_dt] , dtype = float) 我会收到以下错误：
TypeError: only size-1 arrays can be converted to Python scalars
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "...", line 73, in <module>
x = odeint(model, x0, t)
File ".../python3.7/site-packages/scipy/integrate/odepack.py", line 244, in odeint
int(bool(tfirst)))
File "...", line 65, in model
dx_dt = np.array([dq_dt , dk_dt] , dtype = float)
ValueError: setting an array element with a sequence.
那么，有没有办法将dx_dt 转换为浮点数，这样它就不会打扰odeint？
【问题讨论】：
-
你想让 "np.array([dq_dt , dk_dt])" 做什么？ numpy 将其视为对象数组或错误并不奇怪，因为它是异构的。
-
@Redoubts 因为
odeint需要一个函数来返回导数的表达式，在这种情况下dq dt和dx dt。 -
model(x0 , t)产生什么？ -
Numpy 只创建矩形数组。您是否期望
np.array([dq_dt , dk_dt])成为 2xN 矩阵？现在它是 L 形的。
标签： python arrays numpy scipy odeint