Sympy dsolve求解含余弦项二阶ODE时内存溢出崩溃如何解决

• 待求解目标为二阶线性周期系数常微分方程：$y'' + (w/c)^2(\cos(2\pi x/d)+q_0)y=0$，需要同时获取解析解与数值解
• 初始方案调用sympy.dsolve实现符号求解，代码如下：

c = 1
w = 1.5
d = 1
q0 = 2
wc = (w / c)**2

x = symbols('x')
y = Function('y')

equation = Eq(y(x).diff(x, x) + wc * (cos(2 * pi * x / d) + q0) * y(x), 0)
y_x = dsolve(equation)

y_x

• 代码运行时Jupyter直接报错：The kernel has died. It will restart automatically（内核已终止，将自动重启）
• 前置测试：相同调用方式求解2个不含$\cos(kx)$项、附带初始条件的常微分方程时可正常返回结果，初步判断异常与余弦周期项相关，无明确修复方向。

该方程是标准马蒂厄（Mathieu）方程，属于周期系数二阶线性ODE范畴，不存在初等函数构成的闭合解析解。
sympy.dsolve默认会穷举所有内置的初等/特殊函数求解规则做匹配，在找不到初等解的场景下不会快速终止，会持续递归遍历无效求解分支，内存占用持续上涨直到耗尽，最终触发Jupyter内核崩溃。
此前可正常运行的常系数二阶ODE存在固定形式的初等解析解，dsolve可以快速匹配到对应求解规则返回结果，不会触发无效遍历逻辑。

符号解析解（特殊函数形式） #

马蒂厄方程的解析解由第一类、第二类马蒂厄特殊函数线性组合构成，不需要让dsolve做全量规则搜索，直接指定求解hint匹配马蒂厄方程求解逻辑即可，可避免无效内存占用：

from sympy import symbols, Function, Eq, cos, pi, dsolve
c = 1
w = 1.5
d = 1
q0 = 2
wc = (w / c)**2

x = symbols('x')
y = Function('y')
equation = Eq(y(x).diff(x, x) + wc * (cos(2 * pi * x / d) + q0) * y(x), 0)
# 指定匹配马蒂厄方程求解规则，跳过其他无效分支遍历
y_x = dsolve(equation, hint='2nd_linear_mathieu')

如果指定hint后依然存在卡顿，可直接基于sympy内置的mathieuc（偶马蒂厄函数）、mathieus（奇马蒂厄函数）手动构造带两个任意常数的通解即可。

数值解 #

如果不需要符号形式的特殊函数解，直接换用数值积分接口求解即可，运行效率高且不会出现内存溢出问题，推荐使用scipy.integrate.solve_ivp实现：

import numpy as np
from scipy.integrate import solve_ivp

# 固定参数
c = 1
w = 1.5
d = 1
q0 = 2
wc = (w / c)**2

# 将二阶ODE转换为一阶微分方程组
def mathieu_ode(x, y):
    # y[0]为y(x)，y[1]为y'(x)
    return [y[1], -wc * (np.cos(2 * np.pi * x / d) + q0) * y[0]]

# 示例求解配置：求解区间x∈[0, 5]，初始条件y(0)=1、y'(0)=0
x_range = (0, 5)
init_cond = [1, 0]
num_sol = solve_ivp(mathieu_ode, x_range, init_cond, dense_output=True)

# 按需采样计算结果
x_sample = np.linspace(0, 5, 1000)
y_sample = num_sol.sol(x_sample)[0]

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