我正在尝试解决以下一对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 ("nnnSIM 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？
当我通过x0定义进行复制时：
In [41]: x0
Out[41]: array([0, 0])
并在起点进行测试计算：
In [42]: model(x0, 0)
/usr/local/bin/ipython3:29: VisibleDeprecationWarning:
Creating an ndarray from ragged nested sequences (which is a
list-or-tuple of lists-or-tuples-or ndarrays with different
lengths or shapes) is deprecated. If you meant to do this, you
must specify 'dtype=object' when creating the ndarray
Out[42]:
array([0,
array([100. , 99.49028159, 97.9663226 , 95.44365884,
91.94800728, 87.51500386, 82.18984026, 76.02680317,
69.08872084, 61.44632264, 53.17751801, 44.36660217,
35.10339685, 25.48233457, 15.60149599, 5.56161002,
-4.53497306, -14.58532495, -24.48698867, -34.139023 ,
-43.44303157, -52.30416587, -60.63209224, -68.34191273,
-75.35503059, -81.59995152, -87.01301249, -91.53903077,
-95.13186646, -97.75489286, -99.38136988, -99.99471662,
-99.58868039, -98.16740047, -95.74536592, -92.34726785,
-88.00774772, -82.7710442 , -76.69054217, -69.8282285 ,
-62.25406016, -54.04525101, -45.28548466, -36.0640614 ,
-26.47498782, -16.61601846, -6.58765929, 3.5078569 ,
13.56761271, 23.48905528, 33.17104177, 42.51487044,
51.42528687, 59.81145498, 67.58788309, 74.67529543,
81.00144031, 86.50182667, 91.12038155, 94.81002171,
97.53313359, 99.26195678, 99.97886703, 99.67655589,
98.35810522, 96.03695581, 92.73677031, 88.49119201,
83.34350191, 77.34617746, 70.56035759, 63.05521944,
54.90727317, 46.19958194, 37.02091515, 27.46484351,
17.62878515, 7.61301246, -2.48037008, -12.54846682,
-22.48863986, -32.19955543, -41.58221687, -50.54097388,
-58.98449759, -66.8267116 , -73.98766951, -80.39436986,
-85.9815004 , -90.69210386, -94.47815861, -97.30106822,
-99.1320549 , -99.95245291, -99.7538988 , -98.53841671,
-96.31839771, -93.11647348, -88.96528564, -83.90715291])],
dtype=object)
In [43]: _.shape
Out[43]: (2,)
model返回一个包含标量0的2元素对象数组和一个(100，(数组。这不是odeint可以使用的！