【fortran】ニュートン法

　今回は非線形方程式

$$e^x = x^2$$

の解を数値計算で解く。今回はニュートン法で求める。

ニュートン法

　ある関数\(y = f(x)\)を考える。\(f(x)=0\)に対応する適当な近似解\(x_0\)が最初に与えられたものとする。

　次に点\((x_0,f(x_0))\)を通る\(y = f(x)\)の接線方程式

$$y – f(x_0) = f'(x_0)(x-x_0) \tag{1}$$

を考える。式(1)の接線と\(x\)軸の交点を\(x_1\)として新たな近似解とすると、 

$$x_1 = x_0 – \frac{f(x_0)}{f'(x_0)}$$

が得られる。同様に、次は点\((x_1,f(x_1))\)を通る\(y = f(x)\)の接線方程式の\(x\)切片\(x_2\)を求め、近似解を更新する。一般に、近似解\(x_{k+1}\)は\(x_{k}\)から

$$x_{k+1} = x_k – \frac{f(x_k)}{f'(x_k)}　\tag{2}$$

と求められる。式(2)を繰り返し計算し、\(|x_{k+1}-x_{k}|\)が十分小さい値になると\(x_{k+1}\)は解に収束したとみなす。

　 ニュートン法は極めて収束が速いが、初期値や\(f(x)\)によっては解が収束しない場合もある。

サンプルプログラム

　\(f(x) = e^x – x^2\)としている。Newton法では導関数が必要なので、\(g(x) = e^x – 2x\)も定義している。

　使い方はターミナルで

とし、実行して\(x_0\)を代入する。

　例えば\(x_0 = 1\)と入力すると、

と出力される。