如何使用Python Scipy求解特殊二阶常微分方程

给定二阶常微分方程 $y'' + 2y' + 2y = \sin(2x)$，初始条件 $y(0)=0$，$y'(0)=0$，以下是用Scipy求解的具体步骤：

步骤1：将二阶ODE转化为一阶方程组 #

Scipy的solve_ivp仅支持直接求解一阶常微分方程组，因此需要做变量替换：

• 令 $u_1 = y$
• 令 $u_2 = y'$

转化后得到一阶方程组：
$$
\begin{cases}
u_1' = u_2 \
u_2' = -2u_2 - 2u_1 + \sin(2x)
\end{cases}
$$
初始条件对应为 $u_1(0)=0$，$u_2(0)=0$。

步骤2：编写Python代码实现求解 #

导入依赖库 #

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

定义微分方程组 #

def ode_system(x, u):
    u1, u2 = u
    du1_dx = u2
    du2_dx = -2*u2 - 2*u1 + np.sin(2*x)
    return [du1_dx, du2_dx]

设置求解参数并调用求解器 #

# 求解的x范围
x_span = (0, 10)
# 初始条件
u0 = [0, 0]
# 生成用于绘图的x点
x_eval = np.linspace(x_span[0], x_span[1], 1000)

# 调用solve_ivp求解
sol = solve_ivp(ode_system, x_span, u0, t_eval=x_eval)

提取结果并可视化 #

# 提取数值解
y_num = sol.y[0]
x_vals = sol.t

# 计算解析解用于对比
y_analytic = (np.exp(-x_vals)*(2*np.cos(x_vals) + 4*np.sin(x_vals)) - 2*np.cos(2*x_vals) - np.sin(2*x_vals)) / 10

# 绘图对比
plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_num, label='数值解', alpha=0.7)
plt.plot(x_vals, y_analytic, label='解析解', linestyle='--')
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('二阶ODE的数值解与解析解对比')
plt.legend()
plt.grid(True)
plt.show()

代码说明 #

• solve_ivp是Scipy中求解初值问题的核心函数，参数依次为：方程组函数、求解区间、初始条件、指定输出的x点（参数名t_eval为通用命名，实际对应我们的自变量x）。
• 代码中同时计算了解析解，用于验证数值解的正确性，绘图结果会显示两者几乎重合，说明求解结果可靠。

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