pypy.com/

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

c言語で非線形方程式の解法(ニュートン法)

こんにちは。

久しぶりに記事を書いていきます。

今日は、ニュートン法という非線形方程式の解法をやっていきたいと思います。

ニュートン法は、直感的にもわかりやすく、コーディングも比較的楽なので、使いやすいと思います。

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

ニュートン法の考え方

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

f:id:SaidaTaisei:20200719211841p:plain

上の図の青色の線の式のx軸との交点の座標、つまり点Eのx座標を求めることを考えます。

まず、x座標の初期値を設定します。

今回は上図のようにx0を初期値にしたとします。

ここで、x = x_0 における青線の傾きを f'(x_0)、点Aのy座標を f(x_0)だとすると、

点Aにおける接線の方程式は

y = f'(x_0)(x-x_0) + f(x_0)

と書くことができます。

この接線とx軸との交点、すなわち点Bのx座標は、

x_1 = x_0 - f(x_0)/f'(x_0)

同様に、点Cにおける接線を求めて、点Dのx座標を求めると、

x_2 = x_1 - f(x_1)/f'(x_1)

となります。

これを繰り返すことで、点Eのx座標を出していきます。

最終的には、更新する幅が小さくなったところなんかで計算を止めます。

ニュートン法のコード

それでは、c言語を使って、コードを書きます。

今回は以下のように書きました。

#include <stdio.h>

double f(double x);
double _f(double x);

int main(void){
  double x_start = 2.0;
  double err = 100.0;
  double x_pre = x_start;

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

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

double f(double x){
  return 5*x*x*x + 10*x*x - 2*x -5;
}

double _f(double x){
  return 15*x*x + 20*x -2;
}

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

1.275510
0.877004
0.717828
0.689623
0.688756
result:0.688756


今回は以下の式の解を求めています。

y = 5x^3+10x^2-2x-5

グラフにすると、以下の図のようになります。
f:id:SaidaTaisei:20200720141806p:plain

グラフからも、実行結果があっていることが確認できます。

今回は以上で終わります。