C言語で二分法

前回の記事 二分法についての説明

mashiroyuya.hatenablog.com

前回のつづきで今回は二分法のソースコードを書きたいと思います。

二分法の流れは

  • 解を挟む範囲を決めて
  • それを半分に!
  • 半分の範囲のどこに解が入っているかの判定
  • もう一回半分➡︎解がどこにあるかの判定➡︎範囲を縮める
  • 何度もやって範囲が小さくなったら終わり

f:id:mashiroyuya:20160507231847j:plain

この流れでこれをそのままC言語で組んでいきます。

自分が書いたコードでは 収束条件として |x0-xn| < EPS となるようにしています。EPSをいじることで収束条件も変えることができます 。 IMAXは繰り返し回数の条件です。IMAXだけ計算を繰り返すと解が収束しなかったとしてプログラムが終了します。 関数 f(x) を変えたいときはmainの外にある double f(double x)の中を変えてください。

具体的なソースコード

logx = 0 の解を求める二分法のコードです。

#include<stdio.h>
#include<math.h>
#include<stdlib.h> //exit関数のために必要

#define IMAX 1000 //計算回数
#define EPS 1.0e-5 //収束条件

double f(double x);

int main(void){

    double x0,xn,xc;
    int i;

    printf("範囲[x0,xn]\nx0=");
    scanf("%lf",&x0);

    printf("xn=");
    scanf("%lf",&xn);

    xc=(x0+xn)/2;

    if(f(x0)*f(xn)<0.0){    //範囲の中に解があれば計算開始
        for(i=0;i<IMAX;i++){
            if( f(x0)*f(xc)<0.0 ){ //範囲を縮める方向の判定
                xn=xc;
            }else x0=xc;

            xc=(x0+xn)/2;

            if( fabs(x0-xc) < EPS) {
                printf("f(x)=0の解は: x=%lf\n",xc); //収束条件を満たせば解を表示
                exit(1);                           //そしてexitでプログラムを終了
            }
        }printf("解が収束しませんでした\n");
    }printf("[x0,xn]には解が存在しません\n");

    return 0;
}
double f(double x){
    return log(x);
}

少しif文を多用しすぎて見づらい気もします。