非線形方程式の解法 〜 ニュートン・ラフソン法 〜
2次方程式のように解の公式(解析解)がある場合は、それを用いて解を求めることができます。しかし次のような方程式の解を求めたい場合はどうすれば良いでしょう。
$$f(x) = x - \cos(x) = 0\ \ \ \ \ \ \ \ \ \ (1)$$
ここでは、非線形方程式の解法で最もよく用いられているニュートン・ラフソン法を使ってこの問題を解くことにします。ニュートン・ラフソン法の特徴は、ある手続きを用いて一発で解を求めるのではなく、同じ手続きを何回か繰り返すことによって次第に解に近付いていくところにあります。この「同じようなことを繰り返し行う」というのはコンピュータの得意とするところです。では具体的な方法を説明しましょう。
まず$x$として適当な値$x_0$を選びます。次に$(x_0,f(x_0))$における$y=f(x)$の接線を求めます。そして接線が$x$軸と交わる点$(y=0)$の$x$座標の値を$x_1$とします。さらに$(x_1,f(x_1))$における接線を求め、$x$軸との交点を$x_1$とします。このように$x_0,x_1,x_2,\cdots$と続けていけば、次第に解に到達するであろうことは、グラフを書いてみるとわかるでしょう。
さて次にこの繰り返しの方法を定式化しましょう。$(x_k,f(x_k))$における接線の方程式は
$$ y-f(x_k)=f’(x_k)(x-x_k) $$
ですから、これが$x$軸と交わる点$(y=0)$における$x_{k+1}$は
$$x_{k+1}=x_k-\displaystyle{\frac{f(x_k)}{f’(x_k)}} $$
となります。この式を繰り返して実行していけばいいわけです。
では、いつまで繰り返せばよいのでしょうか。普通は、$x_k$と$x_{k+1}$がほとんど同じになったところでこの繰り返しを終了します。つまり、$\epsilon$を小さな数として、
$$\biggl| \displaystyle{\frac{x_{k+1}- x_k}{ x_{k+1}}}\biggr| \leq\epsilon \ \ \ \ \ \ \ \ \ \ \ (2)$$
が成立したところで終了します。これで$x_{k+1}$が解に収束したと判断するわけです。ただし、上の条件を満たしたからといって、$x_{k+1}$が$f(x_{k+1})=0$を満たす厳密解であるとは限りません。あくまでも解の近似値(近似解)です。
それでは、式(1)の方程式を解くプログラムを作ってみましょう。
1 program newton_raphson
2 implicit none
3 double precision x0,eps,f,df
4 integer nitr,maxitr,icond
5
6 external f,df
7
8 data eps / 1.0d-8 /
9
10 write(*,'(a,$)') 'x0 = '
11 read(*,*) x0
12 write(*,'(a,$)') 'maxitr = '
13 read(*,*) maxitr
14
15 call newton( f, df, x0, eps, maxitr, nitr, icond )
16
17 if ( icond .eq. 0 ) then
18 write(*,'(a,i4)') 'number of iterations : ', nitr
19 write(*,'(a,1pe15.8)') 'solutions : ', x0
20 else if ( icond .eq. 1 ) then
21 write(*,'(i4,a)') nitr, ' iterations but not converged'
22 end if
23
24 end program newton_raphson
25
26 subroutine newton( func, dfunc, x0, eps, maxitr, nitr, icond )
27 implicit none
28 double precision x0,x1,eps,func,dfunc
29 integer nitr,maxitr,icond
30
31 icond = 1
32 nitr = 0
33 do while( icond .eq. 1 .and. nitr .lt. maxitr)
34 nitr = nitr + 1
35 x1 = x0 -func(x0) / dfunc(x0)
36 if( abs( (x1 - x0) / x1) .le. eps) icond = 0
37 x0 = x1
38 end do
39
40 end subroutine
41
42 function f(x)
43 double precision f,x
44
45 f = x - dcos(x)
46
47 end function
48
49 function df(x)
50 double precision df,x
51
52 df = 1 + dsin(x)
53
54 end function
上記のプログラムでは、メインプログラム内で関数f と dfをexternal宣言(外部手続きの宣言)することにより、サブルーチンに関数を引数で渡せるようしています。このようにすることにより、1つのプログラム内で様々な関数に対して’newton’のサブルーチンを使用することが可能な、汎用性の高いプログラムになります。
○課題○ 上記のプログラムを、初期値をいろいろ変えて実行しなさい。収束しない場合は、maxitrの値を大きくしてみなさい。また、収束条件である$\epsilon$ (eps) の値を変えて、解への収束状況を調べなさい。
○課題○ 次の方程式の解を求めなさい。
$$x^2 - 2\sin x = 0$$
上式は$x=0$も解にもちますが、$x=0$が解の場合、(2)式の収束の判定条件が使えません。このときは異なる収束条件を用いる必要があります。
○課題○次の方程式の解をニュートン・ラフソン法により求めなさい。maxitr=10、初期値として$x_0=1.2$、 $x_0=2.0$とした場合に、各nitrにおける$x$の値を書き出し、それぞれ初期値に対する振る舞いを調べなさい。
$$x^4 - 6x^2 - 11 = 0$$
上式の解は$x=\sqrt{3+2\sqrt5}=2.733520798347724…$ですが、$x_0=1.2$としたときは$x_k$は振動してしまい、解が求まりません。ニュートン・ラフソン法では、初期値が解の近くではないと、このような振る舞いを示す場合があります。一方で$x_0=2.0$としたときには、有効桁数が反復ごとに2倍になる、2次収束が観測できます。
一般にはニュートン・ラフソン法は収束も速く、非線形方程式を解くのに非常に有力な方法です。またここで用いた繰り返し法で解を求めるテクニックは、色々なところで応用が効きます。しかし方程式が複数の解をもつ場合は、ニュートン・ラフソン法でそれらの解を同時に求めることはできません。初期値によってどの解に収束するかが変わります。また、導関数が解析的に求まっていないと、この解法は使えません(微分を数値的に行う方法はあります)。このような制限があることを十分に理解して利用してください。
非線形方程式の解法 〜 二分法 〜
ニュートン・ラフソン法以外にも非線形方程式の解法はあります。以下では、最も簡単な解法の1つである二分法と呼ばれる方法を紹介します。
非線形方程式$f(x)=0$の解を、以下の手順により二分法で求めることができます。
(1) $f(x_{min})\cdot f(x_{max}) < 0$となるような$ x_{min}$、$ x_{max}$の対を見つける。
(2) $x_0 = (x_{min} + x_{max})/2$として、$f(x_0)$を求める。
(3-1) $f(x_{min})\cdot f(x_0) < 0$の場合：解は$ x_{min}$と$x_0$の間にあるので、$ x_{max}=x_0$とおき、ステップ(2)に戻る。
(3-2) $f(x_{min})\cdot f(x_0) > 0$の場合：解は$x_0$と$ x_{max}$の間にあるので、$ x_{min}=x_0$とおき、ステップ(2)に戻る。
(4) こうして、$| f(x_0) |$が十分に小さくなる($< \epsilon$になる)までステップ(2)と(3)を繰り返すことにより、解が求まる。
○自由課題 (余裕があれば、試してみてください)○ 二分法により、次の方程式の解を求めるプログラムを作成しなさい。また、ニュートン・ラフソン法と比べて収束性 (あるいは収束までにかかる nitrの数)がどのくらい違うか調べなさい。
(a) $f(x)=\sin x - \cos x$
(b) $f(x)=x - 2\ln (x+1)$