使用scipy.integrate.RK45求解微分方程的问题及代码故障排查

嘿，我明白你现在卡在哪了——直接调用RK45后打印出来的是个求解器对象，根本不是你想要的数值数组对吧？这是因为RK45属于显式迭代型的求解器，它不会一次性输出整个区间的结果，得你手动驱动它一步步计算，然后把结果收集起来。我来一步步帮你搞定，先从你测试的非耦合方程入手，再扩展到耦合方程组的情况。

首先得修正你现有代码里的几个小问题：

• RK45要求的函数签名是fun(t, y)，你的f(t,x)符合要求，但要注意参数顺序是时间t在前，状态变量在后。
• 你构造RK45对象时的参数顺序错了！正确的参数顺序是：RK45(fun, t0, y0, t_bound, max_step=inf, rtol=1e-3, atol=1e-6, ...)，你把max_step、rtol的位置搞混了，而且e**-6建议写成1e-6更规范。

下面是修正后的完整代码，能正确收集解并绘图：

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

# 定义微分方程：x' = -x
def f(t, x):
    return -x

# 初始化求解器：t0=0，初始值x0=[1]，求解到t_bound=10
solver = RK45(f, t0=0, y0=[1], t_bound=10, rtol=1e-6, atol=1e-6)

# 用来收集结果的列表
t_values = []
x_values = []

# 循环推进求解，直到到达目标时间t_bound
while solver.status == 'running':
    solver.step()  # 计算下一个时间步的解
    t_values.append(solver.t)  # 记录当前时间
    x_values.append(solver.y[0])  # 记录当前x的值（y是数组，取第一个元素）

# 转成numpy数组方便后续处理
t_values = np.array(t_values)
x_values = np.array(x_values)

# 打印前5个结果示例
print("前5个时间点的解：")
for t, x in zip(t_values[:5], x_values[:5]):
    print(f"t={t:.4f}, x={x:.4f}")

# 绘制解曲线
plt.figure(figsize=(8, 4))
plt.plot(t_values, x_values, label='x(t) = e^(-t)')
plt.xlabel('t')
plt.ylabel('x')
plt.title('非耦合微分方程x\'=-x的解')
plt.legend()
plt.grid(True)
plt.show()

假设你的耦合方程组是：
$x' = f(x, y, t)$
$y' = g(x, y, t)$

这里用一个具体的耦合例子（$x' = -x + y$，$y' = x - 2y$）来演示完整的求解和绘图流程：

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

# 定义耦合微分方程组：状态变量数组[y1, y2]对应x和y
def coupled_fun(t, y):
    x, y_var = y  # 解包状态变量，用y_var避免和函数参数y冲突
    dx_dt = -x + y_var  # x' = -x + y
    dy_dt = x - 2 * y_var  # y' = x - 2y
    return [dx_dt, dy_dt]

# 初始条件：t0=0，x0=1，y0=0
t0 = 0
y0 = [1, 0]
t_bound = 10

# 初始化求解器
solver = RK45(coupled_fun, t0=t0, y0=y0, t_bound=t_bound, rtol=1e-6, atol=1e-6)

# 收集结果
t_values = []
x_values = []
y_values = []

while solver.status == 'running':
    solver.step()
    t_values.append(solver.t)
    x_values.append(solver.y[0])
    y_values.append(solver.y[1])

# 转成numpy数组
t_values = np.array(t_values)
x_values = np.array(x_values)
y_values = np.array(y_values)

# 绘制时间序列解曲线
plt.figure(figsize=(10, 5))
plt.plot(t_values, x_values, label='x(t)')
plt.plot(t_values, y_values, label='y(t)')
plt.xlabel('t')
plt.ylabel('x, y')
plt.title('耦合微分方程组的解')
plt.legend()
plt.grid(True)
plt.show()

# 绘制相图（x vs y），观察变量间的关系
plt.figure(figsize=(6, 6))
plt.plot(x_values, y_values)
plt.xlabel('x')
plt.ylabel('y')
plt.title('耦合方程组的相图')
plt.grid(True)
plt.show()

• RK45是迭代式求解器，必须手动循环调用.step()推进求解，这是它和odeint、solve_ivp（后者可直接返回结果）的核心区别。
• 自定义的微分方程函数必须遵循fun(t, y)的签名，y是状态变量的数组，返回值是对应导数的数组。
• 初始化求解器时，t_bound是你要求解的最终时间，rtol和atol控制求解精度，数值越小精度越高。
• 每次调用.step()后，可通过solver.t获取当前时间，solver.y获取当前状态变量的值。

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