パソコンで数値計算

2分法

数値計算の手法

ある方程式
   (1)
の解の左側と右側にあるx 軸上の点の値をそれぞれa とb とすると必ず次の関係式
   (2)
が成り立ちます。2分法ではこの性質を利用し、方程式の解を求めます。まず上記の式を満たすa; b の初期値として適当に決めておき、そして変数mを
   (3)
とし、これをまず仮の解と見なしてみます。もし初期値として決めたa,bが、正確な解の値からそれぞれ等しい分離れていたら、求めたmの値は正確な解の値となります。もし、
   (4)
なら仮の解は、実際の解より右側にあるということだからbにmの値を代入します。そうでないならば仮の解は、実際の解より左側にあるということなのでaにmの値を 代入します。そして新たに決まったa,bの値を用いてmの値を計算し、また同様に新たなaとbの値を求めていきます。この作業を繰り返していくといずれmの値は正確な値に近づいていくでしょう。そして微小値をepsとしてaとbの値が次の関係式
   (5)
を満たしたら計算終了とし、そこで得られたmを解とします。このように2分法はシンプルな方法で解を求められる方法ですが、解が求まるまでの計算回数が多くあまり効率的ではありません。そこで2分法よりも効率的な方法としてニュートン法と呼ばれる解を求める方法があります。

プログラムソース

今回は方程式 xsin(x)+log(x)=0 を解くとします。プログラム中の関数double f(double x)の内容を変えれば様々な方程式を解くことができます。

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

#define eps 1.0e-5              //収束条件

double f(double x);
void nibun(void); 

int main()
{
        nibun();

        return 0;
}

void nibun(void)
{
        int count;
        double a,b,m;

        count=0;

        printf("範囲の左の値を入力してください.\n");
        scanf("%lf",&a);
        printf("範囲の右の値を入力してください.\n");
        scanf("%lf",&b);

        do {
                count++;
                m=(a+b)/2.0;
                if(f(m)*f(a)<0) b=m;
                else a=m;

//1000回繰り返して解が収束しなかったらプログラムを終了させる
                if(count==1000) {
                        printf("収束しませんでした.\n");
                        exit(1);
                }
        } while (!(fabs(a-b)<eps));             //解が収束したらループを終了

        printf("解の値は %f\n収束するのに %d 回かかりました.",m,count);

}

double f(double x)
{
        return x*sin(x)+log(x);
}

数値計算実行結果

2分法を行う前に、求めたい方程式の解のある程度の予測が必要です。(2)の条件を満たすようにaとbの値をまず初めに決めなければなりません。そのためには、グラフィック系ソフト等を使って方程式の関数を描き、求めた方程式の解のだいたいの範囲を求める必要があります。

範囲の左の値を入力してください.
0.5
範囲の右の値を入力してください.
1.0
解の値は 0.664101
収束するのに 16 回かかりました.

2009/04/19 更新