Heun方法（改进的欧拉法）VS 梯形公式

简短回答：Heun方法（改进的欧拉法）和梯形公式在数学形式上是相同的，但实现方式和数值稳定性不同。

让我详细解释：

---

一、数学形式相同

梯形公式（隐式方法）：

$$y_{n+1}=y_n+\frac{h}{2}\left[f(t_n,y_n)+f(t_{n+1},y_{n+1})\right]$$

这是一个隐式方程，因为右边包含未知的 $y_{n+1}$。

Heun方法（改进的欧拉法）：

1. 预测：${\tilde{y}}_{n+1}=y_n+hf(t_n,y_n)$
2. 校正：$y_{n+1}=y_n+\frac{h}{2}\left[f(t_n,y_n)+f(t_{n+1},{\tilde{y}}_{n+1})\right]$
   把预测代入校正，形式上得到： $$y_{n+1}=y_n+\frac{h}{2}\left[f(t_n,y_n)+f(t_{n+1},y_n+hf(t_n,y_n))\right]$$

看起来和梯形公式一样，但关键区别在于第二个 ( f ) 的参数：

方法	第二个 $f$ 的参数
梯形公式	$f(t_{n+1},y_{n+1})$ ← 未知的真实值
Heun方法	$f(t_{n+1},y_n+hf(t_n,y_n))$ ← 预测值

---

二、本质区别：显式 vs 隐式

1. 梯形公式是隐式的

• 方程：$y_{n+1}=y_n+\frac{h}{2}[f(t_n,y_n)+f(t_{n+1},y_{n+1})]$
• 需要解方程（迭代或牛顿法）
• 对于非线性 $f$，$y_{n+1}$ 无法直接计算

2. Heun方法是显式的

• 先算预测值 ${\tilde{y}}_{n+1}$（显式欧拉）
• 再代入梯形公式的第二个 $f$
• 所有计算都是显式的，可直接计算

---

三、具体例子说明区别

考虑非线性方程：

$$y^{\prime }=y^2,y(0)=0.5,h=0.1$$

用Heun方法（改进欧拉）：

1. 预测：

   $${\tilde{y}}_1=0.5+0.1\times (0.5{)}^2=0.5+0.025=0.525$$

2. 校正：

   $$f_0=(0.5{)}^2=0.25$$ $${\tilde{f}}_1=(0.525{)}^2=0.275625$$ $$y_1=0.5+\frac{0.1}{2}(0.25+0.275625)=0.5+0.05\times 0.525625=0.52628125$$

直接计算，没问题。

用梯形公式（隐式）：

方程：

$$y_1=0.5+0.05[0.25+y_1^2]$$

即：

$$y_1=0.5125+0.05y_1^2$$

或：

$$0.05y_1^2-y_1+0.5125=0$$

这是一个二次方程，需要解：

$$y_1=\frac{1\pm \sqrt{1-4\times 0.05\times 0.5125}}{0.1}=\frac{1\pm \sqrt{0.8975}}{0.1}$$ $$=\frac{1\pm 0.94736}{0.1}=\{\begin{array}{l}19.4736\\ 0.5264\end{array}$$

取物理合理的根：$y_1\approx 0.5264$

需要解方程，得到的结果与Heun方法略有不同：

• Heun: 0.52628125
• 梯形: 0.5264

---

四、数值性质对比

1. 精度阶数

• 两者都是二阶精度（局部截断误差 $O(h^3)$）
• 但常数项不同，所以结果有微小差异

2. 稳定性

这是最重要的区别！

对于测试方程 $y^{\prime }=\lambda y$：

梯形公式（隐式）：

放大因子：

$$R(z)=\frac{1+z/2}{1-z/2},z=h\lambda$$

当 $\text{Re}(\lambda )<0$（衰减问题）时，$|R(z)|<1$ 对所有 $h>0$ 成立。
A-稳定！

Heun方法（改进欧拉）：

放大因子：

$$R(z)=1+z+\frac{z^2}{2}$$

稳定条件：$|1+z+z^2/2|<1$
稳定区域有限，不是A-稳定。

关键结论：梯形公式处理刚性方程更好。

---

五、几何解释

Heun方法（改进欧拉）：

1. 从起点沿切线走到预测点（虚线）

   P: (t_n, y_n) → P': (t_{n+1}, \tilde{y}_{n+1})
       斜率 = f(t_n, y_n)
   

2. 计算终点预测点的斜率

   在P'点：斜率 = f(t_{n+1}, \tilde{y}_{n+1})
   

3. 用起点和预测点斜率的平均

   最终斜率 = (起点斜率 + 预测点斜率)/2
   从起点沿这个平均斜率走
   

梯形公式：

1. 想象一条线，它的平均斜率等于起点斜率和终点斜率的平均
2. 但终点斜率依赖于终点位置
3. 需要调整终点位置使得这个关系成立

就像：

• Heun：先猜终点，用猜的终点斜率
• 梯形：要找真正的终点，使起点和终点的平均斜率正好合适

---

六、从预测-校正角度理解

Heun方法 = PECE的一次迭代

• P（预测）：显式欧拉
• E（求值）：$f(t_{n+1},{\tilde{y}}_{n+1})$
• C（校正）：梯形公式
• E（求值）：$f(t_{n+1},y_{n+1})$（为下一步准备）

只做一次校正，不迭代。

梯形公式的完全迭代

• 预测：某个初始猜测（如上一步值）
• 迭代：反复应用梯形公式 $$y_{n+1}^{(k+1)}=y_n+\frac{h}{2}[f(t_n,y_n)+f(t_{n+1},y_{n+1}^{(k)})]$$
• 直到收敛：$y_{n+1}^{(k+1)}\approx y_{n+1}^{(k)}$

---

七、在实际计算中的表现

简单例子：$y^{\prime }=-10y,y(0)=1,h=0.5$

精确解：$y(0.5)=e^{-5}\approx 0.0067379$

Heun方法：

预测：${\tilde{y}}_1=1+0.5\times (-10)=-4$（已经错了！）
校正：$y_1=1+0.25[(-10)+(-10\times (-4))]=1+0.25[-10+40]=1+7.5=8.5$
完全错误！振荡发散

梯形公式：

方程：$y_1=1+0.25[-10+(-10y_1)]=1-2.5-2.5y_1$
解得：$y_1=-1.5/3.5\approx -0.4286$（至少稳定有界）

原因：$h\lambda =0.5\times (-10)=-5$

• Heun稳定要求：$|1+z+z^2/2|<1$，对 $z=-5$ 不满足
• 梯形公式：A-稳定，总是稳定

---

八、历史与名称

为什么有这么多名字？

1. 梯形公式：来源于数值积分中的梯形法则

   • 把微分方程积分：$y_{n+1}=y_n+{\int }_{t_n}^{t_{n+1}}f(t,y)dt$
   • 用梯形法则近似积分

2. 改进的欧拉法：相对于基本欧拉法的改进

   • 基本欧拉：只用一个斜率
   • 改进欧拉：用两个斜率的平均

3. Heun方法：以德国数学家 Karl Heun 命名

   • 他系统研究了这类预测-校正方法

4. 二阶Adams-Moulton方法：在多步法框架下

   • 这是最简单的隐式多步法

---

九、给初学者的建议

如何区分？

• 如果你直接代入公式计算 → 很可能是Heun方法（改进欧拉）
• 如果你需要解方程 → 是梯形公式（隐式）

实际使用时：

1. 非刚性、简单问题：用Heun方法（计算简单）
2. 刚性、精度要求高：用梯形公式（稳定性好）
3. 现代软件：通常实现为可选项

MATLAB/Python中：

% MATLAB的ode23t使用梯形公式（隐式）
[t, y] = ode23t(fun, tspan, y0);
 
% 自己实现Heun方法
function y_next = heun_step(f, t, y, h)
    k1 = f(t, y);
    y_pred = y + h * k1;
    k2 = f(t + h, y_pred);
    y_next = y + h/2 * (k1 + k2);
end

---

总结

方面	Heun方法（改进欧拉）	梯形公式
数学形式	$y_{n+1}=y_n+\frac{h}{2}[f_n+f(t_{n+1},y_n+h{f}_n)]$	$y_{n+1}=y_n+\frac{h}{2}[f_n+f(t_{n+1},y_{n+1})]$
实现方式	显式：直接计算	隐式：需要解方程
计算量	2次函数求值，无迭代	多次函数求值，需要迭代
稳定性	条件稳定（不是A-稳定）	A-稳定（无条件稳定）
适用问题	非刚性或弱刚性	刚性方程
本质	梯形公式的一次近似	梯形公式的精确实现

关键洞察：

1. 形式上相似，但数值行为不同
2. 显式近似 vs 隐式精确
3. 稳定性是主要区别
4. 刚性方程必须用隐式梯形公式

所以你的观察是对的：它们数学表达式相似，但一个是显式近似，一个是隐式精确解。在数值分析中，这种”形式相似但本质不同”的情况很常见，正是这些微妙区别决定了方法的适用性。