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