使用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/lde/dt-R'（i）/lk-1/lci（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

当你区分一个函数时，你会得到一个函数。所以你需要在某一点上对它求值，才能得到一个数字。要计算symphy表达式，可以使用.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，但希望这种现象是意料之中的，不是问题