关于显式/隐式欧拉及梯形法求解ODE的Python代码正确性问询

问题背景 #

需要数值求解常微分方程：

y′(t)=(0.5−t)(y(t) + 1), y(0) = 1

采用显式欧拉法、隐式欧拉法、梯形法并用Python实现，现有代码结果与同学存在差异，需验证实现正确性。

代码问题分析 #

• 隐式欧拉法实现错误
  隐式欧拉法的核心公式是：
  $$y_{n+1} = y_n + h \cdot f(t_{n+1}, y_{n+1})$$
  对于当前线性ODE，可直接推导得到$y_{n+1}$的解析表达式：
  将$f(t_{n+1}, y_{n+1})=(0.5-t_{n+1})(y_{n+1}+1)$代入公式整理得：
  $$y_{n+1} = \frac{y_n + h(0.5-t_{n+1})}{1 - h(0.5-t_{n+1})}$$
  原代码用显式欧拉的预测值代入$f$计算，并非真正的隐式欧拉，属于近似迭代的错误写法。

• 绘图代码存在未定义变量
  原代码中plt.quiver(T, Y, U, V, angles='xy')一行的T、Y、U、V均未定义，会导致运行报错，需删除该行。

• 显式欧拉法与梯形法实现正确
  显式欧拉的公式$y_{n+1}=y_n + h\cdot f(t_n,y_n)$、梯形法的公式$y_{n+1}=y_n + \frac{h}{2}[f(t_n,y_n)+f(t_{n+1},y_n+h\cdot f(t_n,y_n))]$均正确实现。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 通用求解器
def solveIVP(f, tspan, y0, h, solver):
    t = np.arange(tspan[0], tspan[1]+h, h)
    y = np.zeros(len(t))
    y[0] = y0
    
    for n in range(len(t)-1):
        y[n+1] = solver(f, t[n], y[n], h, t[n+1])  # 传入t_{n+1}给隐式欧拉
    return t, y

# 解析解
def y_function(t):
    return 2 * np.exp(0.5*t - 0.5*t**2) - 1
t_values = np.linspace(-1, 2, 100)
y_analytical = y_function(t_values)

# 显式欧拉法
def explicit_euler(f, tn, yn, h, t_next=None):
    return yn + h * f(tn, yn)

# 隐式欧拉法（修正版）
def implicit_euler(f, tn, yn, h, t_next):
    # 针对当前线性ODE直接求解y_{n+1}
    a = 1 - h*(0.5 - t_next)
    b = yn + h*(0.5 - t_next)
    return b / a

# 梯形法
def trapezoid(f, tn, yn, h, t_next=None):
    k1 = f(tn, yn)
    k2 = f(tn + h, yn + h * k1)
    return yn + h / 2 * (k1 + k2)

# 定义ODE右端函数
def f(t,y):
    return (0.5-t)*(y+1)

# 求解参数
tspan = [0, 1]
y0 = 1
h = 0.2

# 求解各方法
t, y_e = solveIVP(f, tspan, y0, h, explicit_euler)
t, y_i = solveIVP(f, tspan, y0, h, implicit_euler)
t, y_t = solveIVP(f, tspan, y0, h, trapezoid)

# 打印结果表格
print("|   t   | 显式欧拉 | 隐式欧拉 | 梯形法 | 解析解 |")
print("|:-----:|:-------:|:-------:|:-----:|:-----:|")
for n in range(len(t)):
    y_an = y_function(t[n])
    print(f"| {t[n]:0.3f} | {y_e[n]:0.4f} | {y_i[n]:0.4f} | {y_t[n]:0.4f} | {y_an:0.4f} |")

# 绘图
fig, ax = plt.subplots(figsize=(10,6))
plt.plot(t_values, y_analytical, "k-", label="解析解")
plt.plot(t, y_e, "bo-", label="显式欧拉法")
plt.plot(t, y_i, "go-", label="隐式欧拉法")
plt.plot(t, y_t, "ro-", label="梯形法")
plt.xlim(-1, 2)
plt.ylim(-1, 2)
plt.xlabel("$t$", fontsize=16)
plt.ylabel("$y$", fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.legend(fontsize=14)
plt.grid(alpha=0.3)
plt.show()

结果说明 #

修正后，三种数值方法的结果会与解析解更贴近：

• 显式欧拉法为向前一步预测，步长较小时误差可控；
• 隐式欧拉法为向后一步求解，对于该ODE稳定性更好；
• 梯形法为二阶精度方法，结果最接近解析解。

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