如何用scipy.integrate.solve_ivp求解含另一微分方程解的微分方程

你的问题核心是重复求解u的微分方程导致效率极低——每次计算y的导数时都重新跑一遍u的完整求解，完全没必要。正确的做法是先一次性求出u在整个时间区间的解，再通过插值得到任意时刻t对应的u值，代入y的微分方程即可。

具体步骤： #

• 预求解u的微分方程：先一次性算出u在整个时间区间的离散解，保存时间点和对应的u值。
• 构建u的插值函数：用插值方法把离散的u解转换成可以根据任意时刻t实时计算u值的函数，保证数值连续性。
• 代入求解y的微分方程：在y的导数函数中，调用插值函数获取当前t对应的u，再计算导数。

代码实现 #

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d

# 假设已知的常数矩阵/初始条件（根据你的实际情况替换）
A = np.array([[0, 1], [-1, 0]])
B = np.array([[0], [1]])
Omega = np.array([0.1])
S = np.array([[0.5]])
t0, tf = 0, 10
y0 = np.array([1, 0])
u0 = np.array([0.0])

# 1. 预求解u的微分方程
def eqn2(t, u):
    dudt = Omega + u @ S @ u
    return dudt

# 求解u，得到时间点t_u和对应的u值u_sol
u_sol = solve_ivp(eqn2, (t0, tf), u0, method='LSODA', t_eval=np.linspace(t0, tf, 1000))
t_u = u_sol.t
u_vals = u_sol.y.T  # 转置后形状为(n_time_points, u_dim)，方便插值

# 2. 构建u的插值函数
# 如果u是向量/矩阵，对每个维度分别插值
u_interp = interp1d(t_u, u_vals, kind='cubic', axis=0, fill_value="extrapolate")
# kind可选'linear'/'cubic'，根据精度需求选择；fill_value处理超出时间区间的情况

# 3. 定义y的微分方程，使用插值后的u
def eqn1(t, y):
    # 获取当前t对应的u值
    u = u_interp(t)
    # 计算导数（注意矩阵维度匹配，根据你的实际情况调整）
    dydt = (A + B @ u) @ y
    return dydt

# 求解y的微分方程
y_sol = solve_ivp(eqn1, (t0, tf), y0, method='LSODA', t_eval=np.linspace(t0, tf, 1000))

# 查看结果
print("y的解：", y_sol.y)

原方法的问题： #

• 重复计算：每次调用eqn1都重新求解整个u的微分方程，相当于把u的求解跑了上千次（solve_ivp内部会多次调用eqn1），效率极低。
• 数值匹配问题：用np.where(u.t == t)匹配时间点时，由于浮点数精度问题，可能找不到完全相等的t，导致索引错误或警告，插值方法则能避免这个问题。

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