ニュートン・ラプソン法とは？｜非線形方程式の数値解法 – 高校物理からはじめる工学部の物理学

方程式や微分方程式の近似解を数値的に見出す様々な手法がありますが、今回はその中でも代表的なニュートン・ラプソン法と具体例について説明します。

ニュートン・ラプソン法とは？

以下に示す手法で $f(x)=0$ の近似解を数値的に計算する方法をニュートン・ラプソン法と呼ぶ。

$\RM{STEP}1$：適当な初期値 $x=x_0$ を設定する。

$\RM{STEP}2$：以下の式を計算し、近似解 $x_1$ を求める。

\begin{split}
x_1=x_0-\ff{f(x_0)}{f'(x_0)}
\end{split}

$\RM{STEP}3$：$|x_1-x_0|$ が閾値より小さいかどうかを判定する。閾値より小さければ計算終了。閾値より大きい場合は、閾値より小さくなるまでステップ２の計算を繰り返し行い $x_n$ を得る。

$\RM{STEP}4$：$x_n$ を $f(x)=0$ の近似値として出力する。

始めに、ニュートン・ラプソン法の数学的な背景について説明します。

クリックしてジャンプ

ニュートン・ラプソン法とは？
1. ニュートン・ラプソン法における誤差の評価
ニュートンラプソン法を用いたケプラー方程式の数値計算
多変数関数とニュートン・ラプソン法

ニュートン・ラプソン法とは？

冒頭で説明したように、方程式や微分方程式の近似解を数値的に求める手法の一つとして、ニュートン・ラプソン法があります。例えば、解析的には解けないケプラー方程式や、マシュー方程式などの解を近似的に求めたい場合などに使われます。

この節では、ニュートン・ラプソン法の概要とその数学的な背景について説明します。

まず、関数 $y=f(x)$ について考えます。この関数が $x=X$ にて $f(X)=0$ となるとします。この $X$ の位置は下図のようになります。

このとき、任意の点 $x=x_0$ での接線の方程式は、

\begin{split}
y=f(x_0)+f'(x_0)(x-x_0)
\end{split}

と表せます。この接線と $x$ 軸との交点 $x_1$ は、

\begin{split}
x_1=f(x_0)-\ff{f'(x_0)}{f(x_0)}
\end{split}

となります。さらに $x=x_1$ での接線を考え、同様に $x$ 軸との交点 $x_2$ を考えると、

\begin{split}
x_2=f(x_1)-\ff{f'(x_1)}{f(x_1)}
\end{split}

と同じ式が現れます。以下、同様の計算を繰り返すと $f(x)=0$ の解、$X$ を求められることが上図から判断できます。このような手順で近似解を求める手法をニュートン・ラプソン法と呼び、その手順は次のようにまとめられます。

ニュートン・ラプソン法とは？

以下に示す手法で $f(x)=0$ の近似解を数値的に計算する方法をニュートン・ラプソン法と呼ぶ。

$\RM{STEP}1$：適当な初期値 $x=x_0$ を設定する。

$\RM{STEP}2$：以下の式を計算し、近似解 $x_1$ を求める。

\begin{split}
x_1=x_0-\ff{f(x_0)}{f'(x_0)}
\end{split}

$\RM{STEP}4$：$x_n$ を $f(x)=0$ の近似値として出力する。

以下に、ニュートン・ラプソン法により近似解が得られる理由について、その数学的な背景を説明します。

ニュートン・ラプソン法における誤差の評価

ニュートン・ラプソン法における誤差の評価を行います。

まず、$f(x)=0$ の解を $X$ として、$k$ 番目に得られた近似解を $x_k$ とします。このときの誤差 $r_k$ は、

\begin{split}
r_k=x_k-X
\end{split}

と表せます。ここで $f(x_k)$ と $f'(x_k)$ を $x=X$ 周りでテイラー展開してみると、

\begin{split}
f(x_k)&=f(X+r_k) \EE
&=f(X)+f'(X)r_k+\ff{1}{2}f^{”}(X)r_k^2+\cdots \EE
&\NEQ f'(X)r_k+\ff{1}{2}f^{”}(X)r_k^2\,\,\,\,(\because f(X)=0)
\end{split}

\begin{split}
f'(x_k)&=f'(X+r_k) \EE
&=f'(X)+f^{”}(X)r_k+\ff{1}{2}f^{”’}(X)r_k^2+\cdots \EE
&\NEQ f'(X)+f^{”}(X)r_k
\end{split}

次にニュートン・ラプソン法の反復式 $x_{k+1}=x_k-\DL{\ff{f(x_k)}{f'(x_k)}}$ を $r_k$ の式として変形すると、

\begin{split}
x_{k+1}&=x_k-\DL{\ff{f(x_k)}{f'(x_k)}} \EE
(r_{k+1}+X)&=(r_{k}+X)-\DL{\ff{f(x_k)}{f'(x_k)}} \EE
\therefore\,\,r_{k+1}&=r_{k}-\ff{f(x_k)}{f'(x_k)}
\end{split}

となります。これに先程のテイラー展開の結果を適用すると、

\begin{split}
r_{k+1}&=r_{k}-\ff{f'(X)r_k+\ff{1}{2}f^{”}(X)r_k^2}{f'(X)+f^{”}(X)r_k} \EE
&=\ff{1}{2}\ff{f^{”}(X)}{f'(X)+f^{”}(X)r_k}r_k^2
\end{split}

今、近似値が $X$ に十分近ければ、$r_k$ を $0$ と見なせるため、分母を $f'(X)+f^{”}(X)r_k\NEQ f'(X)$ と近似できます。以上より、

\begin{split}
r_{k+1}&=\ff{1}{2}\ff{f^{”}(X)}{f'(X)}r_k^2
\end{split}

とあるステップでの誤差の変化の式が得られました。

上式から分かるように誤差が、$\DL{\left|\ff{1}{2}\ff{f^{”}(x_k)}{f'(x_k)}r_k^2\right|}<1$ のような良い性質を満たすとき、$r_{k+1}$ は $r_{k}$ より小さくなることが分かります。これより、十分な回数反復計算を繰り返せば、近似解が $X$ に収束することが分かります。

ニュートンラプソン法を用いたケプラー方程式の数値計算

早速、ニュートン・ラプソン法を使ってみましょう。今回は、天体力学の重要な式でありながら、超越方程式のために解析的に解けない、ケプラー方程式の近似解を求める方法を考えます。

始めに、ケプラー方程式は、次のような方程式のことです。ただし、$M$ を平均近点離角、$E$ を離心近点離角、$e$ を離心率とします。

\begin{eqnarray}
M = E \,-\, e\sin E
\end{eqnarray}

この方程式の厳密解は、ベッセル関数を用いて

\begin{eqnarray}
E &=& M + 2\sum_{k=1}^{\infty} \ff{1}{k} J_k(ek) \sin kM
\end{eqnarray}

と無限級数の形で表現されます。例えば、$e=0.5, M=60^{\circ}$ として $k=30$ まで計算を実行し、小数点 $4$ 桁までを表示すると次のようになります。

\begin{eqnarray}
E &\NEQ& 88.6398^{\circ}
\end{eqnarray}

この結果をニュートン・ラプソン法により求めた近似解と比較してみます。

さて、ケプラー方程式をニュートン・ラプソン法で解く準備として、方程式を次のように変形します。

\begin{eqnarray}
f(E)=E-e\sin E-M
\end{eqnarray}

こうすることで、$f(E)=0$ となる $E$ は自動的にケプラー方程式を満たすように設定できます。これより、解くべき反復式は、

\begin{split}
E_{k+1}&=E_k-\ff{f(E_k)}{f'(E_k)}\EE
&=E_k-\ff{E_k-e\sin E_k-M}{1-e\cos E_k} \EE
&=\ff{e(\sin E_k-E_k\cos E_k)+M}{1-e\cos E_k}
\end{split}

とできます。それでは、$e=0.5, M=60^{\circ}$ として、$E$ を求めてみましょう。初期値を $E=60^{\circ}$ としてエクセルなどで反復式を実行すると、結果は下表のようになります。

\begin{array}{c|c|c}
ステップ\,k & E_k & 誤差\,|E_k-E_{k+1}| \\
\hline
0 & 60.0000^{\circ} & – \\
\hline
1 & 93.0797^{\circ} & 0.577350 \\
\hline
2 & 88.7235^{\circ} & 7.60300\times 10^{-2} \\
\hline
3 & 88.6399^{\circ} & 1.46060\times 10^{-3} \\
\hline
4 & 88.6398^{\circ} & 5.39600\times 10^{-7} \\
\hline
5 & 88.6398^{\circ} & 7.34968\times 10^{-14} \\
\end{array}

表から分かるように、$4$ ステップ程度で十分な精度の近似解が得られることが分かります。

多変数関数とニュートン・ラプソン法

$\RM{GPS}$ を利用して位置計算を行う場合などでは、多変数関数を解く必要があります。このような多変数関数に対してもニュートン・ラプソン法が適用できると便利と言えます。そこで、多変数関数に対してもニュートン・ラプソン法を適用する方法を考えてみます。

ここでは、以下のような $n$ 次の連立方程式を解くことを考えます。

$$
\left\{
\begin{split}
&f_1(x_1,x_2,\cdots,x_n)=0 \EE
&f_2(x_1,x_2,\cdots,x_n)=0 \EE
&\qquad\qquad\vdots \EE
&f_n(x_1,x_2,\cdots,x_n)=0 \EE
\end{split}
\right.
$$

今、上の連立方程式を全て満足する各変数を $X_1,X_2,\cdots X_n$ として、その点周りのテイラー展開を実行すると、

$$
\left\{
\begin{split}
0&=f_1(X_1,X_2,\cdots,X_n)+\left.\ff{\del f_1}{\del x_1}\right|_{x_1=X_1}(x_1-X_1)\\
&\qquad+\cdots+\left.\ff{\del f_1}{\del x_n}\right|_{x_n=X_n}(x_n-X_n)+O(x^2) \EE
0&=f_2(X_1,X_2,\cdots,X_n)+\left.\ff{\del f_2}{\del x_1}\right|_{x_1=X_1}(x_1-X_1)\\
&\qquad+\cdots+\left.\ff{\del f_2}{\del x_n}\right|_{x_n=X_n}(x_n-X_n)+O(x^2) \EE
&\qquad\qquad\vdots \EE
0&=f_n(X_1,X_2,\cdots,X_n)+\left.\ff{\del f_n}{\del x_1}\right|_{x_1=X_1}(x_1-X_1)\\
&\qquad+\cdots+\left.\ff{\del f_n}{\del x_n}\right|_{x_n=X_n}(x_n-X_n)+O(x^2) \EE
\end{split}
\right.
$$

ここで $f_i(\xi_1,\xi_2,\cdots,\xi_n)=f_i$ と簡略表示をして整理しつつ、行列の表現も用いると

$$
\begin{pmatrix}
f_1 \\
f_2 \\
\vdots \\
f_n \\
\end{pmatrix}
+
\begin{pmatrix}
\left.\ff{\del f_1}{\del x_1}\right|_{x_1=X_1} & \dots & \left.\ff{\del f_1}{\del x_n}\right|_{x_n=X_n} \\
\left.\ff{\del f_2}{\del x_1}\right|_{x_1=X_1} & \dots & \left.\ff{\del f_2}{\del x_n}\right|_{x_n=X_n} \\
\vdots & \ddots & \vdots\\
\left.\ff{\del f_n}{\del x_1}\right|_{x_1=X_1} & \dots & \left.\ff{\del f_n}{\del x_n}\right|_{x_n=X_n} \\
\end{pmatrix}
\begin{pmatrix}
x_1-X_1 \\
x_2-X_2 \\
\vdots \\
x_n-X_n \\
\end{pmatrix}
\begin{split}
=\B{0}
\end{split}
$$