使用SciPy求解ODE时无法调用已定义函数运算的问题

我是一名大一数学与物理专业学生，几乎没有Python使用经验，这是我首次使用SciPy处理常微分方程（ODE）。编写代码时遇到了问题：无法在另一个ODE函数dxdt中调用之前定义的dvdt函数进行除法运算。我尝试将dvdt的返回代码直接替换到dxdt中，但因v在dxdt内未定义而失败；尝试将dvdt转为float类型也不行，因为无法将函数转为浮点型。

原代码 #

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
import math

def dvdt(t, v):
    
    k = 0.3
    Fs = 200000
    m0 = 10000
    lam = 200
    
    return (Fs + k*v**2) / (m0 - lam*t)

v0 = 0

def dxdt(t, x):
    
    k = 0.3
    Fs = 200000
    m0 = 10000
    lam = 200
    
    # 问题出在这里
    return math.sqrt(Fs / (dvdt*k))

# 把dvdt的返回代码直接替换进来也没用，因为dxdt里没定义v :/

x0 = 0

t = np.linspace(0, 30, 100)

sol_v = integrate.solve_ivp(dvdt, t_span = (0, max(t)), y0=[v0], t_eval=t)
v_sol = sol_v.y[0]

sol_x = integrate.solve_ivp(dxdt, t_span = (0, max(t)), y0=[x0], t_eval=t)
x_sol = sol_x.y[0]

plt.plot(t, v_sol)
plt.ylabel('$v(t)$', fontsize=22)
plt.xlabel('$t$', fontsize=22)

plt.plot(t, x_sol)
plt.ylabel('$x(t)$', fontsize=22)
plt.xlabel('$t$', fontsize=22)

plt.show()

问题分析 #

原代码中dxdt函数里直接写dvdt*k是错误的，因为dvdt是一个函数对象，不能直接和数值相乘。而直接替换dvdt的返回表达式又会因为v在dxdt作用域内未定义而报错。

修正后代码 #

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
import math

def dvdt(t, v):
    
    k = 0.3
    Fs = 200000
    m0 = 10000
    lam = 200
    
    return (Fs + k*v**2) / (m0 - lam*t)

v0 = 0

def dxdt(t, x):
    
    k = 0.3
    Fs = 200000
    m0 = 10000
    lam = 200
    
    # 正确调用dvdt函数，传入当前t和v=0的参数
    return math.sqrt(Fs / (dvdt(t, 0)*k))

x0 = 0

t = np.linspace(0, 30, 100)

sol_v = integrate.solve_ivp(dvdt, t_span = (0, max(t)), y0=[v0], t_eval=t)
v_sol = sol_v.y[0]

sol_x = integrate.solve_ivp(dxdt, t_span = (0, max(t)), y0=[x0], t_eval=t)
x_sol = sol_x.y[0]

plt.plot(t, v_sol, label="v(t)")
plt.plot(t, x_sol, label="x(t)")

plt.ylabel('速度', fontsize=15)
plt.xlabel('时间', fontsize=15)

plt.legend(["v(t) - 速度", "x(t) - 高度"])
plt.title("火箭运动状态")
plt.show()

修正说明 #

修正后的代码在dxdt中正确调用了dvdt函数，传入了当前的时间t和初始速度0作为参数，这样就能得到dvdt在该时刻的数值结果，进而完成后续的计算。同时对图表的标签和标题进行了中文本地化，更易读。

内容的提问来源于stack exchange，提问作者Raiden