中點法的圖示說明,假設
等於
的實際值。中點法計算
讓紅色的弦近似平行於中點處的切線(綠色)
數值分析裡的中點法(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.