【fortran】ニュートン法で方程式の解を求めてみた!

プログラミング

今回は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)\)の解を求める方法です。

詳しく言うと

  1. 初期値\(x_0\)を決め、その点の傾き\(f'(x_0)\)を求める。
  2. 傾きの直線 (灰色の線) と\(x\)軸(黒の横線)との交点\(x_1\)を求める。
  3. 手順1の\(x_0\)を\(x_1\)に変えて1.2を繰り返す。
  4. 繰り返していき、\(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)\)が0に近づくまでループを繰り返します。

 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
が出ることが確認できます。

タイトルとURLをコピーしました