使用 scicpy.integrate.odeint 和 sympy 符号时出错

Question

我正在尝试解决以下系统：d²i/dt² + R'(i)/L di/dt + 1/LC i(t) = 1/L dE/dt 作为一组耦合的一阶微分方程：

• di/dt = k
• dk/dt = 1/L dE/dt - R'(i)/L k - 1/LC i(t)

这是我正在使用的代码：

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

#Define model: x = [i , k] 

def RLC(x , t):
    
    i = sp.Symbol('i')
    t = sp.Symbol('t')

    #Data: 
    E = sp.ln(t + 1)
    dE_dt = E.diff(t)
    
    R1 = 1000   #1 kOhm
    R2 = 100    #100 Ohm  
    R = R1 * i + R2 * i**3
    dR_di = R.diff(i)
    
    i = x[0]
    k = x[1]
    L = 10e-3   #10 mHy
    C = 1.56e-6 #1.56 uF
    
    #Model
    di_dt = k
    dk_dt = 1/L * dE_dt - dR_di/L * k - 1/(L*C) * i
    dx_dt = np.array([di_dt , dk_dt])
    
    return dx_dt

#init cond:
x0 = np.array([0 , 0])

#time points:
time = np.linspace(0, 30, 1000)

#solve ODE:
x = odeint(RLC, x0, time)

i = x[: , 0]

但是，我收到以下错误：TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

所以，我不知道sympy 和odeint 是否不能很好地协同工作。或者可能是因为我将t 定义为sp.Symbol 造成的问题？

Answer 1

当你区分一个函数时，你会得到一个函数。因此，您需要在某个时间点对其进行评估才能获得一个数字。要评估一个 sympy 表达式，您可以使用 .subs()，但我更喜欢 .replace()，它感觉更强大（至少对我而言）。

您必须尝试让每个变量都有自己的名称以避免混淆。例如，您从一开始就将浮点输入 t 替换为 sympy Symbol，从而丢失了 t 的值。变量 x 和 i 也在外部范围内重复，如果它们的含义不同，这不是一个好的做法。

以下内容应避免混淆，并有望产生您所期望的结果：

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


# Define model: x = [i , k]

def RLC(x, t):
    # define constants first
    i = x[0]
    k = x[1]
    L = 10e-3  # 10 mHy
    C = 1.56e-6  # 1.56 uF
    R1 = 1000  # 1 kOhm
    R2 = 100  # 100 Ohm

    # define symbols (used to find derivatives)
    i_symbol = sp.Symbol('i')
    t_symbol = sp.Symbol('t')

    # Data (differentiate and evaluate)
    E = sp.ln(t_symbol + 1)
    dE_dt = E.diff(t_symbol).replace(t_symbol, t)

    R = R1 * i_symbol + R2 * i_symbol ** 3
    dR_di = R.diff(i_symbol).replace(i_symbol, i)
    
    # nothing should contain symbols from here onwards
    # variables can however contain sympy expressions

    # Model (convert sympy expressions to floats)
    di_dt = float(k)
    dk_dt = float(1 / L * dE_dt - dR_di / L * k - 1 / (L * C) * i)
    dx_dt = np.array([di_dt, dk_dt])

    return dx_dt


# init cond:
x0 = np.array([0, 0])

# time points:
time = np.linspace(0, 30, 1000)

# solve ODE:
solution = odeint(RLC, x0, time)

result = solution[:, 0]
print(result)

需要注意的是：i = x[0] 的值似乎在每次迭代中都非常接近于 0。这意味着dR_di 基本上一直停留在1000。我不熟悉 odeint 或您的特定 ODE，但希望这种现象是意料之中的，不是问题。