Newton-Raphson法を用いた円周率近似 - それが僕には楽しかったんです。

Newton-Paphson法の基礎原理

Newton-Raphson法とは次の式を満たすxの値を見つけるアルゴリズム
${ \displaystyle f(x) = g(x) - C = 0 }$
ここではg(x)を任意の微分可能な関数としてCを実数定数とする。
つまりNewton-Rapshon法は次の式を満たすxの値を見つけるアルゴリズムとなる。

${ \displaystyle g(x) = C }$

つまり以下の様にg,Cが与えられた時は2の平方根を求めることになる。

${ \displaystyle g(x) = x^2, C = 2 }$

実際の計算

基礎原理はうえで紹介したので実際に解を近似して求める。というのも上記の式を満たすxがすぐに正確に見つかるわけではないことによる。
そこで与えられたf(x)を一回微分して接線の方程式を求め、その接線とx-軸の接点におけるx座標で再びf(x)に戻り微分し接点のx座標を求める…。
これを再帰的に繰り返すことで求めたいxに近似していくこととなる。
その漸化式は次のようになる。

${ \displaystyle x_{n+1} = x_n - \frac{f(x_n)}{\frac{df}{dx}(x_n)} }$

この漸化式をC言語を用いたコードとして表現していく。
ここでは平方根を求めることにする。

コード化する

ひとまず次のプログラムをベースにする。

#include<stdio.h>

int main(){
    return 0;
}

ここで関数fとそれを一回微分したdfを関数として作る。

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

//プロトタイプ宣言
double Func(double x,double C);//変数xとCが必要となる
double dFunc(double x);//fを微分しているのでxのみ必要

int main(){
    return 0;
}

double Func(double x, double C){
    return x*x - C; //x^2-Cを計算して返す
}

double dFunc(double x){
    return 2.0*x; //2xを返す
}

これで準備は出来た。

次に微分して…。となるがCで微分を扱うのは難しいので数学的に考えてfとxのみの式に直す。
それは以下になる。

${ \displaystyle \frac{df}{dx}(x_n) = \frac{f(x_n)}{x_n-x_{n+1}} }$

というわけでNewton-Raphson法を用いて平方根を求める。ここでは誤差が ${10e-6}$ になるまで続けることにする。

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

//追加した
#define EPS 10e-6

//プロトタイプ宣言
double Func(double x,double C);//変数xとCが必要となる
double dFunc(double x);//fを微分しているのでxのみ必要

double Newton(double C);//平方根を求めたい値を実数で受け取る

int main(){

    printf("%12.10f\n",Newton(2) );

    return 0;
}

//Newton-Rapshon法を用いてcの平方根を求める。
double Newton(double C){
    
    double x = C;
    double dx = -Func(x,C)/dFunc(x);
    
    while( fabs(dx) > EPS){
        x += dx;
        printf("%12.10f\n",x);
        dx = -Func(x,C)/dFunc(x);
    }
    return x;
}

double Func(double x, double C){
    return x*x - C; //x^2-Cを計算して返す
}

double dFunc(double x){
    return 2.0*x; //2xを返す
}

これで平方根をNewton-Raphson法を使って求めることができる。
ただし、誤差があることに注意。

円周率近似

円周率近似も漸化式を用いて行う。
http://oto-suu.seesaa.net/article/291513400.html
上記のサイトがそれを紹介しているが半径が ${\frac{1}{2}}$ なので、半径を単位円において考えることにする。
基本的には正n角形の周を求め円周に近似していくことになる。

その場合、漸化式は次のようになる。
${ \displaystyle L_{n+1} = \sqrt{ 2-2\sqrt{1-(\frac{L_n}{2})^2} } or L_{n+1} = \sqrt{2-\sqrt{4-(L_n)^2}} }$

ソースコードにする

これをソースコードにするとこうなる

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

#define EPS 1.0e-8
#define POW 2

double Newton(double);
double Func(double, double);
double dFunc(double);

int main(){

    double Ln = Newton(2.0);//L_1 = sqrt(2)
    printf("%12.10f\n",Ln);
    int n = 1;
    int pow_base = 2;
    int n_max = 20;

    while( n < n_max ){ 
        
        printf( "%12.10f\n", pow_base * Ln );
        Ln = Newton( 2-Newton(4 - pow(Ln, 2)) );
        n += 1;
        pow_base *= POW;

    }

    return 0;
}


//平方根を求める
double Newton(double C){
    
    double x = C;
    double dx = -Func(x, C) / dFunc(x);
    
    while( fabs(dx) > EPS ){
        x += dx;
        dx = -Func(x, C) / dFunc(x);
    }

    return x;
}


double Func(double x, double C){
    return x*x-C;
}

double dFunc(double x){
    return 2.0 * x;
}

これで円周率の近似値が求められた。