今回はfortranを使い、ニュートン法で方程式の解を求めたいと思います。
このようなコードになります。
program newton implicit none external f, fd integer :: i real(8) :: x0, x real(8) f, fd ! !--------- data --------------------------- write(*,*) ' initial x = ' read(*,*) x0 !--------- do ----------------------------- do i=1,100 x=x0-f(x0)/fd(x0); if(abs(f(x))<=1.d-8) then write(*,*) 'root=',x exit end if x0=x end do end program newton !-------- function ------------------------- function f(x) implicit none real(8) :: x,f f = x**3-10*x**2+23*x-15 end function f function fd(x) implicit none real(8) :: x,fd fd = 3*x**2-20*x+23 end function fd
ニュートン法で方程式の解を求める
ニュートン法とは
ニュートン法とは上の図のようなイメージで\(f(x)\)の解を求める方法です。
詳しく言うと
- 初期値\(x_0\)を決め、その点の傾き\(f'(x_0)\)を求める。
- 傾きの直線 （灰色の線） と\(x\)軸（黒の横線）との交点\(x_1\)を求める。
- 手順1の\(x_0\)を\(x_1\)に変えて1.2を繰り返す。
- 繰り返していき、\(f(x)\)（赤色の線）の値が0に十分近くなったら、その\(x\)は方程式の解とみなせる。
といった感じです。
手順2の傾きの直線（灰色の直線）の式は
$$
y-f(x_0) = f'(x_0)(x-x_0)
$$
となります。
\(y=0\)のとき、\(x=x_1\)となるので、代入して変形すると
$$
x_1 = x_0 – \frac{f(x_0)}{f'(x_0)}\tag{1}
$$
となり\(x=x_1\)が求まります。
- この\(x_1\)を\(f(x)\)に代入し、\(f(x_1)\)が十分0に近かければ、\(x_1\)が答えになります。
- もし違えば、この\(x_1\)を(1)式に代入し、\(f(x)\)が0に近くなるまでこの操作を繰り返します。
コード解説
ではコードの解説をしていきます。
今回のコードでは
$$
f(x) = x^3-10x^2+23x-15
$$
の解を求めています。
初期値
9~10行目では初期値\(x_0\)を10に設定しています。
x0=10 write(*,*) ' initial x = ',x0
繰り返し処理
12~19行目では先ほどの繰り返し処理が行われています。
13行目で(1)式を使って、傾きの直線と\(x\)軸との交点\(x_1\)を求めています。
x=x0-f(x0)/fd(x0);
そして14~17行目のif文で「\(x\)を代入したときに、\(f(x)\)が\(1.0\times10^{-8}\)よりも小さければ（0に十分近ければ）繰り返しを終える(exit)」という命令を行っています。
if(abs(f(x))<=1.d-8) then write(*,*) 'root=',x exit end if
この条件が満たされなければ、\(x_0\)に\(x\)が代入されて（18行目）\(f(x)\)が０に近づくまでループを繰り返します。
x0=x
関数
23~33行目は解を求めたい関数
$$
f(x) = x^3-10x^2+23x-15
$$
とその導関数
$$
f'(x) = 3x^2-20x+23
$$
を書いています。
答え
このコードを動かしてみると
initial x = 10.000000000000000 root= 7.0329344278299191となります。
これで解の一つが求めることができました！！
補足
本当は解が三つあるはずですが、一つしか出ないことに疑問がある方もいると思います。
どういうことかというと、今回は初期値\(x_0=10\)に最も近い解が出ています。
例えばこれを\(x_0=1\)などとするとまた異なる答え
1.2226743196073671が出ることが確認できます。
