pypy.com/

Python、Unity、FX自動化などを勉強しています。あと、コーラと車も好きです。そこらへんについて、たまに記事を書きます。

c言語で非線形方程式の解法（セカント法）

数値計算

こんにちは。

今回は、セカント法という非線形方程式の解法について書いていきたいと思います。

セカント法は、以前に書いたニュートン法と考え方が似ているので、以前書いたニュートン法の記事も参考になると思います。

以下にリンクを貼っておきますので、気になる方は是非ご覧ください。

taiboman.com

それではやっていきたいと思います。




セカント法の考え方

以下の図を使って、説明していきたいと思います。

上図の青線のx軸との交点を求めることを考えます。

まず、セカント法では、初期値を2つ設定します。

今回は、上図のように x_1、x_2 を初期値として設定したとします。

青線の点Bにおける接線の傾きを考えると、点Aと点Bの座標から以下のように近似できます。

f'(x_2)≃\frac{f(x_2)-f(x_1)}{x_2-x_1}

これを用いて、点Cのx座標 x_3 を求めると、

x_3=x_2-\frac{f(x_2)}{f'(x_2)}

となります。

同様にして、

f'(x_3)≃\frac{f(x_3)-f(x_2)}{x_3-x_2}

であり、

x_4=x_3-\frac{f(x_3)}{f'(x_3)}

となります。

これを繰り返すことで、青線とx軸の交点を求めることができ、それが、方程式の解となります。

セカント法のc言語コード

今回は以下の式の解を求めることを考えます。

f(x)=sin(x)+3cos(x)-2

グラフだと下図のものになります。



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

ニュートン法の時とほとんど同じですね。

収束の判定は更新幅が一定値を下回ったらということにしました。

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

double f(double x);
double _f(double x, double x_pre);

int main(void){
  double x_start = 2.0;
  double err = 100.0;
  double x_pre1 = x_start+0.01;
  double x_pre2 = x_pre1;

  while(err>0.000001){
      x_start = x_start - f(x_start)/_f(x_start,x_pre2);
      err = (x_start - x_pre1) * (x_start - x_pre1);
      x_pre2 = x_pre1;
      x_pre1 = x_start;
      printf("%lf\n",x_start);
  }

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

double f(double x){
  return sin(x)+3.0*cos(x)-2;
}

double _f(double x, double x_pre){
  return (f(x)-f(x_pre)) / (x-x_pre);
}

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

1.255593
1.215637
1.207973
1.207828
result:1.207828

上で示したグラフの交点とほぼ一致していることがわかります。

今回は以上です。