2 非線型方程式

2.1 概要

非線形方程式²

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

の解の値が必要になることが，工学の問題でしばしばある．工学の問題では，数学でやったよう な厳密な解の必要はなく，精度の良い近似解で良いことが多い．近似解といっても， $10^{-10}$程度の精度のことを言っており，この程度の近似解が必要となる．

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

この授業では4通りの計算テクニックを学習したが，重要³なのは

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

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)のみで ある．計算に必要な式は分かったが，数列がどのように真の解$\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 double func(double x);
   5 double dfunc(double x);
   6 
   7 /*================================================================*/
   8 /*       main function                                            */
   9 /*================================================================*/
  10 int main(){
  11   double eps=1e-15;         /* precision of calculation */
  12   double x[IMAX+10];
  13   char temp;
  14   int i=-1;
  15 
  16   printf("\ninitial value x0 = ");
  17   scanf("%lf%c", &x[0], &temp);
  18 	
  19   do{
  20     i++;
  21     x[i+1]=x[i]-func(x[i])/dfunc(x[i]);
  22 
  23     printf("  %d\t%e\n", i, x[i+1]);
  24 
  25     if(fabs((x[i+1]-x[i])/x[i])<eps) break;
  26   }while(i<=IMAX);
  27 
  28   if(i>=IMAX){
  29     printf("\n  not converged !!! \n\n");
  30   }else{
  31     printf("\niteration = %d  solution  x = %20.15f\n\n",i,x[i+1]);
  32   }
  33 
  34   return(0);
  35 }
  36 
  37 /*================================================================*/
  38 /*       define function                                          */
  39 /*================================================================*/
  40 double func(double x){
  41   double y;
  42 
  43   y=x*x*x-3*x*x+9*x-8;
  44 
  45   return(y);
  46 }
  47 
  48 /*================================================================*/
  49 /*       define derived function                                  */
  50 /*================================================================*/
  51 double dfunc(double x){
  52   double dydx;
  53 
  54   dydx=3*x*x-6*x+9;
  55 
  56   return(dydx);
  57 }

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
著者: 山本昌志