如何用Python SciPy求解含离散时变系数的复值耦合微分方程？

核心问题分析 #

你遇到的问题源于ODE求解器（无论是odeintw还是SciPy的complex_ode）采用自适应步长算法，求解过程中会查询任意时间点t的系数值，而非仅你指定的100个离散点。当前代码用int(t)//freq索引离散数组，仅能处理整数倍freq的t，对非整数点会出错，必须把离散的H矩阵转换成连续可查询的函数。

方案1：用插值实现连续的H系数函数 #

利用NumPy的插值函数，把离散的(100,2,2)数组m转换成可接受任意t输入的函数。线性插值计算快，适合多数时变系统；若系统变化剧烈，可改用三次样条插值提升精度。

代码示例：

import numpy as np
from odeintw import odeintw

# 加载离散H矩阵数组m，形状(100,2,2)
m = np.load('your_h_matrix.npy')

dt = 0.1
freq = 100
num_time_steps = 100
# 生成与m索引一一对应的离散时间点数组
t_discrete = np.array([step * freq * dt for step in range(num_time_steps)])

# 定义插值函数，对H矩阵的四个元素分别做线性插值
def get_H(t):
    H11 = np.interp(t, t_discrete, m[:,0,0])
    H12 = np.interp(t, t_discrete, m[:,0,1])
    H21 = np.interp(t, t_discrete, m[:,1,0])
    H22 = np.interp(t, t_discrete, m[:,1,1])
    return H11, H12, H21, H22

def equations(z, t):
    z1, z2 = z
    H11, H12, H21, H22 = get_H(t)
    dz1dt = z1 * H11 + z2 * H12
    dz2dt = z1 * H21 + z2 * H22
    return [dz1dt, dz2dt]

# 显式定义复数初始条件，避免类型错误
z0 = np.array([0+0j, 1+0j])
# 求解目标时间点
t_eval = t_discrete

# 求解ODE
solution = odeintw(equations, z0, t_eval)

z1, z2 = solution[:, 0], solution[:, 1]

方案2：使用SciPy的complex_ode求解 #

complex_ode原生支持复值耦合方程，只需注意导数函数的参数顺序（t在前，z在后），同样搭配插值处理H系数：

代码示例：

import numpy as np
from scipy.integrate import complex_ode

# 加载m数组并生成对应离散时间点
m = np.load('your_h_matrix.npy')
dt = 0.1
freq = 100
num_time_steps = 100
t_discrete = np.array([step * freq * dt for step in range(num_time_steps)])

def get_H(t):
    H11 = np.interp(t, t_discrete, m[:,0,0])
    H12 = np.interp(t, t_discrete, m[:,0,1])
    H21 = np.interp(t, t_discrete, m[:,1,0])
    H22 = np.interp(t, t_discrete, m[:,1,1])
    return H11, H12, H21, H22

def equations(t, z):
    z1, z2 = z
    H11, H12, H21, H22 = get_H(t)
    dz1dt = z1 * H11 + z2 * H12
    dz2dt = z1 * H21 + z2 * H22
    return [dz1dt, dz2dt]

# 初始化求解器
solver = complex_ode(equations)
solver.set_initial_value(np.array([0+0j, 1+0j]), t_discrete[0])

# 逐步求解并保存结果
solution = np.zeros((num_time_steps, 2), dtype=np.complex128)
solution[0] = solver.y
for i in range(1, num_time_steps):
    solver.integrate(t_discrete[i])
    solution[i] = solver.y

z1, z2 = solution[:, 0], solution[:, 1]

关键注意事项 #

• 插值方式选择：若H矩阵变化平缓，np.interp足够；若系统有剧烈突变，可改用scipy.interpolate.CubicSpline做三次样条插值。
• 复数类型声明：初始条件需显式定义为复数（如0+0j），避免求解器自动转换为实数类型导致计算错误。
• 时间点对应：确保t_discrete与m的索引严格一一对应，插值时才不会出现错位。

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