如何用数组中的时变参数求解微分方程组（ODEs）

算力阿峰

2026-7-19

问题：在odeint中随时间步传入时变系数求解微分方程组

我用odeint求解了含ne、Zeff、Epar常数系数的二阶微分方程组，现在想改成每个时间步从数组中读取对应系数值来求解，但尝试结合args参数失败了，原代码如下：

import numpy as np
from scipy.integrate import odeint

def f(Q, t):
    # constants
    ne=5E19
    ln=17
    permit=8.8542E-12
    m_e=9.1094E-31
    c=2.9979E8
    e=1.6E-19
    B0=6
    R0=0.935
    Zeff=3
    Epar=0.2
    Er=(e**3*ne*ln)/(4*np.pi*permit**2*m_e*c**2)
    D = Epar/Er
    alpha = 1+Zeff
    Fgy = (2*permit*B0**2)/(3*ne*m_e*ln)
    Fgc = Fgy*((m_e*c)/(e*B0*R0))**2
    nur=ne*e**4*ln/(4*np.pi*permit**2*m_e**2*c**3)
    

    # assign each ODE to a vector element
    qp, q = Q
    gamma=(1+q**2)**(1/2)
    qpp2=q**2-qp**2
    # define each ODE
    dqpdt = (D-gamma*(gamma+alpha)*qp/(q**3)-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3*qp/q)*nur
    dqdt = (D*qp/q-gamma**2/q**2-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3)*nur

    return [dqpdt,dqdt]

解决方案

odeint本身不支持自动传入时变参数，但可以通过时间插值的方式，在每个求解步根据当前时间获取对应的系数值，具体实现如下：

1. 准备时变系数数组

先定义和求解时间数组t长度一致的时变系数数组，例如：

# 示例时间数组（假设求解0到1e-6秒，共1000步）
t = np.linspace(0, 1e-6, 1000)
# 示例时变系数数组（模拟随时间变化的ne、Zeff、Epar）
ne_arr = 5E19 + 1E19 * np.sin(t * 1e7)
Zeff_arr = 3 + np.cos(t * 1e7)
Epar_arr = 0.2 + 0.1 * np.sin(t * 8e6)

2. 修改微分方程函数

在函数中添加时变系数数组和时间数组作为参数，用np.interp根据当前时间t插值得到对应时刻的系数：

def f(Q, t, t_grid, ne_arr, Zeff_arr, Epar_arr):
    # 固定常数（不变的参数保留在这里）
    ln=17
    permit=8.8542E-12
    m_e=9.1094E-31
    c=2.9979E8
    e=1.6E-19
    B0=6
    R0=0.935

    # 根据当前时间t，从数组中插值得到对应系数
    ne = np.interp(t, t_grid, ne_arr)
    Zeff = np.interp(t, t_grid, Zeff_arr)
    Epar = np.interp(t, t_grid, Epar_arr)

    # 计算依赖于时变系数的中间变量
    Er=(e**3*ne*ln)/(4*np.pi*permit**2*m_e*c**2)
    D = Epar/Er
    alpha = 1+Zeff
    Fgy = (2*permit*B0**2)/(3*ne*m_e*ln)
    Fgc = Fgy*((m_e*c)/(e*B0*R0))**2
    nur=ne*e**4*ln/(4*np.pi*permit**2*m_e**2*c**3)
    

    # assign each ODE to a vector element
    qp, q = Q
    gamma=(1+q**2)**(1/2)
    qpp2=q**2-qp**2
    # define each ODE
    dqpdt = (D-gamma*(gamma+alpha)*qp/(q**3)-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3*qp/q)*nur
    dqdt = (D*qp/q-gamma**2/q**2-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3)*nur

    return [dqpdt,dqdt]

3. 调用odeint时传入额外参数

使用args参数把时变系数相关的数组传入函数：

# 初始条件
Q0 = [0.1, 0.5]  # qp初始值，q初始值
# 调用odeint
sol = odeint(f, Q0, t, args=(t, ne_arr, Zeff_arr, Epar_arr))

# sol是求解结果，每行对应一个时间步的qp和q值
qp_sol = sol[:, 0]
q_sol = sol[:, 1]

关键注意事项

确保时间数组t_grid（即传入的t）是严格单调递增的，否则np.interp无法正确工作
如果你的系数是阶梯式突变（不是连续变化），可以用np.searchsorted替代np.interp，直接定位到当前时间对应的数组索引：
```
idx = np.searchsorted(t_grid, t, side='right') - 1
ne = ne_arr[idx]
Zeff = Zeff_arr[idx]
Epar = Epar_arr[idx]
```

内容的提问来源于stack exchange，提问作者Eduardo Gordo Quiroga