非線形方程式（二分法）

/*********************************/
/* 二分法による exp(x)-3x=0 の根 */
/*      coded by Y.Suganuma      */
/*********************************/
#include <stdio.h>
#include <math.h>

double bisection(double (*)(double), double, double, double,
                 double, int, int *);
double snx(double);

int main()
{
	double eps1, eps2, x, x1, x2;
	int max, ind;
/*
          データの設定
*/
	eps1 = 1.0e-10;
	eps2 = 1.0e-10;
	max  = 100;
	x1   = 0.0;
	x2   = 1.0;
/*
          実行と結果
*/
	x = bisection(snx, x1, x2, eps1, eps2, max, &ind);

	printf("   ind=%d  x=%f  f= %f\n", ind, x, snx(x));

	return 0;
}

/****************/
/* 関数値の計算 */
/****************/
double snx(double x)
{
	double y;
	y = exp(x) - 3.0 * x;
	return y;
}

/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解                  */
/*      fn : f(x)を計算する関数名                        */
/*      x1,x2 : 初期値                                   */
/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
/*      eps2 : 終了条件２（｜f(x(k))｜＜eps2）           */
/*      max : 最大試行回数                               */
/*      ind : 実際の試行回数                             */
/*            （負の時は解を得ることができなかった）     */
/*      return : 解                                      */
/*********************************************************/
#include <math.h>

double bisection(double(*f)(double), double x1, double x2,
                 double eps1, double eps2, int max, int *ind)
{
	double f0, f1, f2, x0 = 0.0;
	int sw;

	f1 = (*f)(x1);
	f2 = (*f)(x2);

	if (f1*f2 > 0.0)
		*ind = -1;

	else {
		*ind = 0;
		if (f1*f2 == 0.0)
			x0 = (f1 == 0.0) ? x1 : x2;
		else {
			sw = 0;
			while (sw == 0 && *ind >= 0) {
				sw    = 1;
				*ind += 1;
				x0    = 0.5 * (x1 + x2);
				f0    = (*f)(x0);

				if (fabs(f0) > eps2) {
					if (*ind <= max) {
						if (fabs(x1-x2) > eps1 && fabs(x1-x2) > eps1*fabs(x2)) {
							sw = 0;
							if (f0*f1 < 0.0) {
								x2 = x0;
								f2 = f0;
							}
							else {
								x1 = x0;
								f1 = f0;
							}
						}
					}
					else
						*ind = -1;
				}
			}
		}
	}

	return x0;
}