中点法的图示说明,假设
等于
的实际值。中点法计算
让红色的弦近似平行于中点处的切线(绿色)
数值分析里的中点法(midpoint method)是求解常微分方程的一种数值方法,属于单步法。
![{\displaystyle y'(t)=f(t,y(t)),\quad y(t_{0})=y_{0}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7c7aaa582de7b3049f85f8b50ee6785dc70c7996)
上式的显式中点法为
![{\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f(t_{n},y_{n})\right),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5bbf9b1957cd660a5a41bb8b0a6e7f647f690bb1) | | 1e |
隐式中点法为
![{\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},{\frac {1}{2}}(y_{n}+y_{n+1})\right),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c285ea6dccef176c1b8da3eb4a635ec5d8c2736a) | | 1i |
。此处
是步阶长度,是一个小的正数,
,而
是计算
的近似值。显式中点法也称为是改良的欧拉方法(modified Euler method)[1],隐式中点式是最简单的配置法(Collocation method),可用在哈密顿力学,也就是辛积分器。
此方法的名称是因为上述公式会用到函数在
位置的值,也就是在
以及
中点时的值,前者的值已知,后者的值未知,在计算中点的值时,也需要有一些假设,才能进行计算。
配合图示(见右图)会比较容易理解此一方法。在原始的欧拉方法里,会用
,计算曲线在
处的切线。此切线和垂直线
的交点即为下一个点
的值。不过,若此函数的二阶导数在时间
和
的区间均为正, 或是均为负(如图中的例子),随着
加大,曲线和切线的距离会越来越远,导致大误差。图中所画的中点的切线(上方,绿线)比较可以近似这段曲线。不过因为要求解的就是此一曲线,时间
的资讯已知,其他资讯不足,可能无法精准的画出中点处的切线。
因此,会先用原始的欧拉方法估计在中点的
值,再用
以及估计的中点资讯,计算切线斜率。用改善后的切线从
计算
的值。最后一步即为图的红色弦线。因为红色弦线是估计值,用的是
在中点处的估计值,不一定真的会和绿线(真正的中点切线)平行,仍会有误差。
中点法每一步的局部误差是
,全域误差是
。其运算比欧拉方法要大,但在
的过程中,中点法的误差会比欧拉方法降低的更快。
此法也是高阶方法(如龙格-库塔法)的范例之一。
微分方程
的数值积分求解图。 蓝:欧拉法
红:精确解:
![{\displaystyle y=e^{t}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d4383514b086d37fd2316fd8baec32739b094291)
步长
。
同一个方程的数值方法求解,步长缩减为
,可以看出中点法的收敛比欧拉法要快
中点法可以视为是改良版的欧拉方法
![{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}),\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f25a25532bf6b3cc9ef75c48f0156e1527a3182f)
因此用类似的方式推导。
推导欧拉法的关键是以下的近似等式
![{\displaystyle y(t+h)\approx y(t)+hf(t,y(t))}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dcf16939499cc3e20c9109427a0c0b56ccf269fd) | | 2 |
是源自以下的斜率公式
![{\displaystyle y'(t)\approx {\frac {y(t+h)-y(t)}{h}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/72bd3e22562f069752d5741ea6db4a1debbac352) | | 3 |
,其中
在中点法中,将(3)改为更准确的式子
![{\displaystyle y'\left(t+{\frac {h}{2}}\right)\approx {\frac {y(t+h)-y(t)}{h}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8011d5154d073c85421545d5f73e65f40fff8b0b)
因此可以得到以下类似(3)式,计算
的式子
![{\displaystyle y(t+h)\approx y(t)+hf\left(t+{\frac {h}{2}},y\left(t+{\frac {h}{2}}\right)\right).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0897744a7a89d22bdc557f7544462e5f8cc252c7) | | 4 |
在上式中,因为还不知道
在
的值,无法用上式直接计算
。解法是用欧拉方法来求解
:
![{\displaystyle y\left(t+{\frac {h}{2}}\right)\approx y(t)+{\frac {h}{2}}y'(t)=y(t)+{\frac {h}{2}}f(t,y(t)),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/78bc7903ae0d3e9e2f3bddfe45002b00123bdfb9)
代入(4)式中,可得
![{\displaystyle y(t+h)\approx y(t)+hf\left(t+{\frac {h}{2}},y(t)+{\frac {h}{2}}f(t,y(t))\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2a65c85ee6e02a762f8e7ad86fab7dd953521519)
以上则是显式的中点法(1e)。
隐式中点法(1i)可以将此时间内的波形假设为直线,中间步数
的值用
和
的平均值来表示
![{\displaystyle y\left(t+{\frac {h}{2}}\right)\approx {\frac {1}{2}}{\bigl (}y(t)+y(t+h){\bigr )}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/eeb47c2967490a8ba029aae5de700f0da51adb6d)
因此
![{\displaystyle {\frac {y(t+h)-y(t)}{h}}\approx y'\left(t+{\frac {h}{2}}\right)\approx k=f\left(t+{\frac {h}{2}},{\frac {1}{2}}{\bigl (}y(t)+y(t+h){\bigr )}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fba04b08b10ebd565754afb108e49fde3c0fa247)
因为隐式方法的时间对称性,局部误差中
的偶次项都消去了,局部误差的阶数是\
。
- Griffiths, D. V.; Smith, I. M. Numerical methods for engineers: a programming approach. Boca Raton: CRC Press. 1991: 218. ISBN 0-8493-8610-1.
- Süli, Endre; Mayers, David, An Introduction to Numerical Analysis, Cambridge University Press, 2003, ISBN 0-521-00794-1 .
- Burden, Richard; Faires, John. Numerical Analysis. Richard Stratton. 2010: 286. ISBN 978-0-538-73351-9.