ニュートン・ラフソン法 [fortran]

こんにちは、さくらえびです。
本記事ではニュートン・ラフソン法(Newton-Raphson Method)について、説明とプログラムを記載します。
ニュートン・ラフソン法は、方程式の求根アルゴリズムの一種です。
求根アルゴリズムとは、厳密定義はさておき、簡単に説明すると、\(f(x)=0\)となる\(x\)を求めるアルゴリズムです。
ニュートン・ラフソン法の式の導出方法は別の本で確認してください。
式は下記の数列が成り立ちます。
\[
x_{n+1} = x_(n) – \frac{f(x)}{f'(x)}
\]
そして、\(f(x)=0\)が成り立つとき、\(x_{n}\)が収束するというものです。
まずはプログラム。微分関数を用いて、こんな風に作ります。

subroutine NEW_LAP(fun,x_ini,ans)
implicit none
external fun
double precision DIF_F,fun,x_ini,ans
double precision xnew,xold
xnew = x_ini
do while (abs(xnew - xold) >1.0E-8)
xold = xnew
xnew = xold - fun(xold)/DIF_F(fun,xold)
enddo
ans = xnew
end subroutine

本プログラムの仕様では、\(10^{-8}\)まで収束すれば良いとします。
ではわかりやすく、下記の式について、試してみましょう。
\[
f(x) = x^2-2
\]
\(f(x) = 0\)が成り立つときは、直感的に、\(x = \pm{\sqrt{2}} \)とわかるため、
\(x = 1.41421356…\)という値になってほしい。
では、計算してみましょう。
テストコードは下記のようにします。

program test
implicit none
external func
double precision x,func,ans,err
call NEW_LAP(func,5.d0,ans)
err = ans-dsqrt(2.d0)
print *,ans,dsqrt(2.d0),err
end program
function func(x)
implicit none
double precision func,x
func = x**2-2
end function

結果は下記のとおり。

x<0側から始めると\(- \sqrt{2} \)に収束していることがわかります。