この非線形方程式は、図1のようにのx軸と交わる点 に実数解を持つ。ここだけとは限らないが、少なくともこの交わる点は解である。この点 の値は、コンピューターを用いた反復(ループ)計算により探すことができる。
この授業では4通りの計算テクニックを学習したが、重要3なのは
実際の数値計算は、 であるような2点 から出発する。 そして、区間を2分する点に対して、を計算を行う。 ならばをと置き換え、 ならばをと置き 換える。絶えず、区間の間に解があるようにするのである。この操作 を繰り返して、区間の幅が与えられた値 よりも小さく なったならば、計算を終了する。解へ収束は収束率1/2の一次収束である。
実際にこの方法で
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 }
まずは、この数列の漸化式を求める。関数上の点 の接 線を引き、それとx軸と交点 である。まずは、 を求める ことにする。点 を通り、傾きが の直線の方 程式は、
(4) |
計算の終了は、
(6) |
ニュートン法を使う上で必要な式は、式(5)のみで ある。計算に必要な式は分かったが、数列がどのように真の解に収束 するか考える。 と真値の差の絶対値、ようするに誤差を 計算する。 を忘れないで、テイラー展開を用いて、計算を進める と
ニュートン法の特徴をまとめると次のようになる。
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 }
二分法に比べて、ニュートン法は収束が早く良さそうであるが、次に示すように解へ収束 しない場合があり問題を含んでいる。
非線形方程式
(8) |
このようにニュートン法は解に収束しないで、振動する場合がある。こうなる と、プログラムは無限ループに入り、永遠に計算し続ける。これは資源の無駄 遣いなので、慎むべきである。通常は、反復回数の上限を決めて、それを防ぐ。 ニュートン法を使う場合は、この反復回数の上限は必須である。
実際には収束しない場合のほうが稀であるので、ニュートン法は非常に強力な 非線型方程式の解法である。ただ、反復回数の上限を決めることを忘れてはならない。