パソコンで数値計算

ニュートン法

数値計算の手法

方程式の解を数値計算で求めるためにニュートン法と呼ばれる方法があり、この方法では接線の性質を利用します。ある方程式f(x)=0があるとします。y = f(x) として関数の曲線をx,y平面に書きます。描かれた曲線とx軸が交わる点が解であり、方程式によって解の個数が1つであったり、複数あったりします。

そこで1つの解に注目し、近辺のx軸上の点を適当に選びその値をa とします。そしてx=aにおける曲線上の点で接する接線を描き、その接線とx軸上の交点はある場合を除いて実際の解に近づくことがわかると思います。その交点のx軸上の値をnewaとすると、f(newa)=0より
   (1)
となり、この式からnewaについて求めると
   (2)
となります。このように求めたnewaの値をaに代入して、同様に(2)から次のnewaの値を求めます。するとどんどんnewaの値は正しい解の値に近づいていくでしょう。そこで正しい解が求まったとする条件式を
   (3)
としてこの式を満たした場合、計算を終わりにしnewaの値を解とします。この条件式は、正しい解の近くにaの値をとったほど、newaの値との差が小さくなることを利用しています。epsは微小な値とし、この値を小さくするほど正確な解の値が求まります。

このニュートン法は2分法と比べると、だいぶ速く解が求まります。しかし、曲線の形によっては上手く解が求まらない場合もありますので注意が必要です。例えば、もし接線がx軸と水平になってしまう、すなわち接線の傾きの値が0になってしまった場合、(2)式の右辺は発散してしまい解は求められません。また曲線が複雑な形をしていると接線の傾きがあちらこちら変化して、上手く解が求まらない事もあります。よって数値計算を行って出てきた値が実際に解が正しいかどうか描いた曲線から判断することも必要です。

プログラムソース

ここではxsin(x)+log(x)=0 の解を求めす。左辺の関数を描きある程度の解の値の予測が必要です。そうしないと計算結果の解の値が正しい解の値を出しているかどうかわからないからです。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define max 1000                //最大繰り返し回数
#define eps 1.0e-5              //収束条件
 
double f(double x);
double df(double x);
void newton(void);
 
int main()
{
        newton();
 
        return 0;
}
 
void newton(void)
{
        int count;
        double a,newa;
 
        count=0;
 
        printf("初期値を入力してください。\n");
        scanf("%lf",&a);
 
        for(;;) {
                count++;
 
                newa=a-f(a)/df(a);
 
                if(fabs(newa-a)<eps) break;
 
                a=newa;
 
                if(count==max) {
                        printf("収束しませんでした。\n");
                        exit(1);
                }
        }
 
        printf("解の値は %f\n収束するのに %d 回かかりました。\n",
            newa,count);
 
}
 
double f(double x)
{
        return x*sin(x)+log(x);
}
 
double df(double x)
{
        return sin(x)+x*cos(x)+1.0/x;
}


数値計算実行結果

プログラムを実行すると、初期値の入力を求められますので、解に近い値を入力します。

初期値を入力してください。
0.5
解の値は 0.664102
収束するのに 4 回かかりました。

2009/04/19 更新