2012年2月23日木曜日

Newton-Roaphson法による方程式の求解

関数のグラフ化により数値解の個数とおおよその解が判別出来たところで、
その解を求めていく作業に入ります。 

全部書くのは大変なのでwikipediaさんのお力を借ります。
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

ニュートン法はwikipediaにもあるように、 まず、任意のx座標(wikiではxn)を決めて、
そこのy値での接線を引きます。



接線を引くと接線の傾きが0でなければ、x軸のどこかと交差します。(wikiではxn+1の点) 
そのxn+1の点を再びxnとして用い、そこのy値の接線を引きます。
そしてまた接線との交点をxn+1とします。 これの繰り返しで方程式の近似解が得られます。
x軸と接線との角をθとすると、tanθは

   

ここで右の二式で求まっていないものはxn+1だけなのでそれを変形すると、

 

となるので現在分かっているxn,f(xn),f'(xn)から未知の点xn+1を求められます。
 これを繰り返すと設定した点xnから大体一番近い解に収束していきます。

 Newton-Roaphson法を使うにあたって必要なこと。
 1.方程式が判っている。(関数のグラフ化では0.5X^3 + 0.75X^2 - 2X - 2)
 2.解の個数が判っている。(同じく3個)
 3.接線の方程式(方程式の一階微分)が判っている。(同じく1.5X^2 + 1.5X - 2)
 4.解の大体の値が判っている。
 こんな感じです。
使用時はValue of X にそこそこ解に近い適当な値(xn)を入れてください。
forの中で10回ループしてますけど、大体そのくらいで収束します。

二つ以上の解の場合、たとえば上の図だと初期値が3以下であれば解は1くらいの値に、
3以上であれば5くらいの値に収束します。
3の場合はエラーだと思います。(接線の傾きが0のためx軸と交わらない)


 以下リスト。

#include<stdio.h>
#include<stdlib.h>

float func(float x);
float dfunc(float x);

int main()
{
 int i;
 float x_init, p_new, p_old;

 printf("Value of X ->");
 scanf("%f", &x_init);

 p_old = x_init;

 for (i = 0; i < 10; i++) {
  p_new = p_old - func(p_old) / dfunc(p_old);
  printf("[%2d] ---- %f\n", i + 1, p_new);
  p_old = p_new;
 }
}

float func(float x)
{
 float value;

 value = 0.5 * x * x * x + 0.75 * x * x - 2.0 * x - 2.0;

 return value;
}

float dfunc(float x)
{
 float value;

 value = 1.5 * x * x + 1.5 * x - 2.0;

 return value;
}

0 件のコメント:

コメントを投稿