2 非線型方程式

2.1 概要

非線形方程式²

$$f(x)=0$$ (1)

の近似解$x$を求める--ことが，ここでの問題である．

この非線形方程式は，図1のように$y= f(x)$のx軸と交わる点 に実数解を持つ．ここだけとは限らないが，少なくともこの交わる点は解である．この点 の値は，コンピューターを用いた反復(ループ)計算により探すことができる．ここの講義 では，次の2通りの方法で近似解を計算した．

1. 2分法
2. ニュートン-ラフソン法(ニュートン法)

図 1: $f(x)=x^3-3x^2+9x-8$の関数．x軸との交点が解である．

2.2 二分法(bisection method)

2.2.1 計算方法

閉区間$[a,\,b]$で連続な関数$f(x)$の値が，

$$f(a)f(b)<0$$ (2)

ならば， $f(\alpha)=0$となる$\alpha$が区間$[a,\,b]$にある．

実際の数値計算は， $f(a)f(b)<0$であるような2点 $a,\,b(a<b)$から出発する． そして，区間$[a,\,b]$を2分する点$c=(a+b)/2$に対して，$f(c)$を計算を行う． $f(c)f(a)<0$ならば$b$を$c$と置き換え， $f(c)f(a)>0$ならば$a$を$c$と置き 換える．絶えず，区間$[a,b]$の間に解があるようにするのである．この操作 を繰り返して，区間の幅$\vert b-a\vert$が与えられた値 $\varepsilon$ よりも小さく なったならば，計算を終了する．解へ収束は収束率$1/2$の一次収束である．

実際にこの方法で

$$x^3-3x^2+9x-8=0$$ (3)

を計算した結果を図2に示す．この図より，$f(a)$と$f(b)$の関係 の式(2)を満たす区間$[a, b]$が1/2ずつ縮小していく様子がわかる．この 方法の長所と短所は，以下の通りである．

長所

閉区間$[a,b]$に解があれば，必ず解に収束する．間違いなく解 を探すので，ロバスト(robust：強靭な)な解法と言われている． 次に示すニュートン法とは異なり，連続であればどんな形の関数 でも解に収束するので信頼性が高い方法と言える．さらに，解の 精度も分かり便利である．解の誤差は，区間の幅$\vert b-a\vert$以下で ある．

短所

収束が遅い(図6)．一次収束である．

図 2: $f(x)=x^3-3x^2+9x-8$の実数解を2分法で解散し，その解の収束の 様子を示している．初期値は $a=-1,\,b=11$として，最初の解$c=x_0=5$が求 まり，順次より精度の良い $x_1,\,x_2,\,x_3,\cdots$が求まる．それが，解析解 $x=1.1659\cdots$(x軸との交点)に収束していく様子が分かる．

2.2.2 アルゴリズム

関数はあらかじめ，プログラム中に書くものとする．更に，計算を打ち切る条 件もプログラム中に書くものとする．そうすると，図 3のような2分法のフローチャートが考えられる．

図 3: 2分法のフローチャート

2.2.3 プログラム

このプログラムをしっかり理解せよ．テストで出題された場合，全く同じプログラムを書 く必要はない．同じアルゴリズムでも同一のプログラムになるとは限らない．

   1 #include <stdio.h>
   2 #define EPS (1.0e-15)  /* precision of calculation */
   3 
   4 double func(double x);
   5 
   6 /*=============================================================*/
   7 /*       main function                                         */
   8 /*=============================================================*/
   9 int main(void){
  10   double a, b, c, test;
  11   char temp;
  12   int i=0;
  13 
  14   do{
  15     printf("\ninitial value a = ");
  16     scanf("%lf%c", &a, &temp);
  17 
  18     printf("initial value b = ");
  19     scanf("%lf%c", &b, &temp);
  20 
  21     test=func(a)*func(b);
  22 
  23     if(test >= 0){
  24       printf("   bad initial value !!  f(a)*f(b)>0\n\n");
  25     }
  26   }while(test >= 0);
  27 
  28   if(b-a < 0){
  29     c=a;
  30     a=b;
  31     b=c;
  32   }
  33 
  34   while(b-a > EPS){
  35 
  36     c=(a+b)/2;
  37     if(func(c)*func(a) < 0){
  38       b=c;
  39     }else{
  40       a=c;
  41     }
  42 
  43     i++;
  44   }
  45 
  46   printf("\nsolution x = %20.15f\n\n",c);
  47 
  48   return 0;
  49 }
  50 
  51 /*=============================================================*/
  52 /*       definition function                                   */
  53 /*=============================================================*/
  54 double func(double x){
  55   double y;
  56 
  57   y=x*x*x-3*x*x+9*x-8;
  58 
  59   return y;
  60 }

2.3 実数解のニュートン法(Newton's method)

2.3.1 計算方法

関数$f(x)$のゼロ点$\alpha$に近い近似値$x_0$から出発する．そして，関数 $f(x)$上の点 $(x_0,f(x_0))$での接線が，x軸と交わる点を次の近似解$x_1$ とする．そして，次の接線がx軸と交わる点を次の近似解$x_2$とする．同じこ とを繰り返して $x_3,\,x_4,\cdots$を求める(図4)．こ の計算結果の数列 $(x_0,\,x_2,\,x_3,\,x_4,\cdots)$は初期値$x_0$が適当で あれば，真の解$\alpha$に収束する．

まずは，この数列の漸化式を求める．関数$f(x)$上の点 $(x_i,\,f(x_i))$の接 線を引き，それとx軸と交点$x_{i+1}$ である．まずは，$x_{i+1}$ を求める ことにする．点 $(x_i,\,f(x_i))$を通り，傾きが $f^\prime(x_i)$の直線の方 程式は，

$$y-f(x_i)=f^\prime(x_i)(x-x_i)$$ (4)

である．$y=0$の時の$x$の値が$x_{i+1}$にである．$x_{i+1}$は，

$$x_{i+1}=x_i-\frac{f(x_i)}{f^\prime(x_i)}$$ (5)

となる．$x_i$から$x_{i+1}$計算できる．これをニュートン法の漸化式と言う． この漸化式を用いれば，次々と近似解を求めることができる．

計算の終了は，

$$\left\vert\frac{x_{i+1}-x_i}{x_i}\right\vert<\varepsilon$$ (6)

の条件を満たした場合とするのが一般的である． $\varepsilon$は計算精度を 決める定数で，非常に小さな値である．これ以外にも計算の終了を決めること は可能である．必要に応じて，決めればよい．実際に式 (3)を計算した結果を図4に 示す．接線との交点が解に近づく様子がわかるであろう．

ニュートン法を使う上で必要な式は，式(5)のみで ある．計算に必要な式は分かったが，数列$x_i$がどのように真の解$\alpha$に収束 するか考える．$x_{i+1}$ と真値$\alpha$の差の絶対値，ようするに誤差を 計算する． $f(\alpha)=0$を忘れないで，テイラー展開を用いて，計算を進める と

$$\begin{aligned}\vert\alpha-x_{i+1}\vert &=\left\vert\alpha-x_i+... ...=\left\vert O\left((\alpha-x_i)^2\right)\right\vert \end{aligned}$$

となる．$i+1$番目の近似値は，$i$番目に比べて2乗で精度が良くなるのである．これを， 二次収束と言い，非常に早く解に収束する．例えば，$10^{-3}$ の精度で近似解が得られ ている場合，式(5)を再び計算するだけで，$10^{-6}$程度の 精度で近似解が得られるということである．一次収束の2分法よりも，収束が早いことが 分かる．

ニュートン法の特徴をまとめると次のようになる．

長所

初期値が適当ならば，収束が非常に早い(図 6)．

短所

初期値が悪いと，収束しない(図7)．収束 しない場合があるので，反復回数の上限を決めておく必要がある．

図 4: $f(x)=x^3-3x^2+9x-8$の実数解をニュートン法で計算し，解の収 束の様子を示している．初期値$x_0=5$から始まり，接線とx軸の交点からよ り精度の高い回を求めている．

2.3.2 アルゴリズム

2分法同様，関数と計算を打ち切る条件はプログラム中に書くものとする．そ うすると，図5のようなニュートン法のフローチャー トが考えられる．

図 5: ニュートン法のフローチャート

2.3.3 プログラム

このプログラムをしっかり理解せよ．テストで出題された場合，全く同じプログラムを書 く必要はない．同じアルゴリズムでも同一のプログラムになるとは限らない．

ニュートン法では，導関数の計算が必要である．通常の方程式であれば，導関数は計算で きるはずである．計算した導関数をプログラム中に記述する．

   1 #include <stdio.h>
   2 #include <math.h>
   3 #define IMAX 50
   4 #define EPS (1.0e-15)            /* precision of calculation */
   5 
   6 double func(double x);
   7 double dfunc(double x);
   8 
   9 /*================================================================*/
  10 /*       main function                                            */
  11 /*================================================================*/
  12 int main(void)
  13 {
  14   double x[IMAX+10];
  15   char temp;
  16   int i=-1;
  17 
  18   printf("\ninitial value x0 = ");
  19   scanf("%lf%c", &x[0], &temp);
  20 	
  21   do{
  22     i++;
  23     x[i+1]=x[i]-func(x[i])/dfunc(x[i]);
  24   }while(i<=IMAX && fabs((x[i+1]-x[i])/x[i]) >= EPS);
  25 
  26   if(i>=IMAX){
  27     printf("\n  not converged !!! \n\n");
  28   }else{
  29     printf("\niteration = %d\nsolution  x = %20.15f\n\n",i,x[i+1]);
  30   }
  31 
  32   return 0;
  33 }
  34 /*================================================================*/
  35 /*       definition function                                      */
  36 /*================================================================*/
  37 double func(double x)
  38 {
  39   double y;
  40 
  41   y=x*x*x-3*x*x+9*x-8;
  42 
  43   return y;
  44 }
  45 /*================================================================*/
  46 /*       definition derived function                              */
  47 /*================================================================*/
  48 double dfunc(double x)
  49 {
  50   double dydx;
  51 
  52   dydx=3*x*x-6*x+9;
  53 
  54   return dydx;
  55 }

2.4 ニュートン法と2分法の比較

2.4.1 解への収束速度

図6に，二分法とニュートン法の解への近 づき具合を示す．二分法に比べ，ニュートン法が解への収束が早いことがわかる．前者は 二次収束で，後者は一次収束であることがグラフより分かる．二分法は，10回の計算で， $2^{-10}=1/1024$程度になっている．

二分法に比べて，ニュートン法は収束が早く良さそうであるが，次に示すように解へ収束 しない場合があり問題を含んでいる．

図 6: 二分法とニュートン法の計算回数(反復回数)と誤差の関係

2.5 ニュートン法の問題点

アルゴリズムから，2分法は解に必ず収束する．ただし，この方法は，収束のスピードが 遅く，それが欠点となっている．一方，ニュートン法は解に収束するとは限らない．初期 条件に依存する場合がある．厳密にその条件を求めるのは大変なので，初期条件により収 束しない実例を示す．

非線形方程式

$$3\tan^{-1}(x-1)+\frac{x}{4}=0$$ (8)

の解を計算することを考える．これは，初期値のより，収束しない場合がある．例えば初 期値$x_0=3$の場合，図7のように収束しない．これを初期値 $x_0=2.5$にすると図8のように収束する．

このようにニュートン法は解に収束しないで，振動する場合がある．こうなる と，プログラムは無限ループに入り，永遠に計算し続ける．通常は反復回数の上限を決め ，無限ループを防ぐ．ニュートン法を使う場合は，この反復回数の上限は必須である．

実際には収束しない場合のほうが稀であるので，ニュートン法は非常に強力な 非線型方程式の解法である．ただ，反復回数の上限を決めることを忘れてはならない．

図 7: ニュートン法で解が求まらない場合

図 8: ニュートン法で解が求まる場合

ホームページ: Yamamoto's laboratory
著者: 山本昌志