使用 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() 感觉更强大（至少对我而言）。

您必须尝试让每个变量都有自己的名称以避免混淆。例如，您从一开始就用 sympy Symbol 替换了 float 输入 t，从而丢失了 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，但希望这种现象是预料之中的并且不是问题。

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

Error using scicpy.integrate.odeint and sympy symbols

python

numpy

scipy

sympy

odeint