c言語で非線形方程式の解法（はさみうち法）

こんにちは。

今回は、はさみうち法という非線形方程式の解法をやっていきたいと思います。

はさみうち法は、以前紹介した二分法を改良したもので、二分法よりも基本的に収束が速い方法となります。

以前書いた二分法の記事は以下より見ることができるので、ぜひご覧ください。

それでは、書いていきます。

はさみうち法の考え方

以下の図を用いて考えます。

はさみうち法では、まず初期値を二つ設定します。

今回は を初期値として設定したとして、考えていきます。

初期値から、 を求めることを考えます。

は点A と点B を通る直線とx軸の交点のx座標であるから、

となります。

ここで、 の正負を判定します。

正であれば、 を用いて を求め、

負であれば、 を用いて を求めることになります。

今回の図の場合だと、 は負であるから、 を求めた時と同様にして、

となります。

これを繰り返すことで、解を求めることができます。

収束の判定は の絶対値が小さくなるxになった時点とします。

はさみうち法のc言語コード

今回は前回の二分法で用いたのと同じ以下の式の解を求めます。

グラフだと、以下のようになります。

コードは以下のように書きました。

#include <stdio.h>
#include <math.h>

double f(double x);

int main(void){
  double x_start1 = 2.0;
  double x_start2 = -1.0;
  double err = 100.0;
  double a = x_start1;
  double b = x_start2;
  double c;

  while(err>0.00000001){
      c = (a*f(b)-b*f(a)) / (f(b)-f(a));
      double fc = f(c);
      err = fc*fc;
      if(fc < 0){
        b = c;
      }
      else{
        a = c;
      }
      printf("%lf\n",c);
  }

  printf("result:%lf\n",c);
}

double f(double x){
  return exp(x)*sin(x)-cos(x);
}

コードを見るとわかるように、二分法のものとほぼ同じになります。

実行結果は以下のようになります。

-0.680697
-0.323821
0.005526
0.249499
0.396097
0.470941
0.505400
0.520417
0.526795
0.529472
0.530591
0.531058
0.531252
0.531333
0.531367
result:0.531367

上のグラフを見る限り、うまく解が求まっているように見えます。

今回は以上です。

他にも数値計算の記事は書いているので、興味のある方はご覧ください。