非線形方程式（セカント法）

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

double secant(double, double, int, double, double, int *,
              double(*)(double));
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 = secant(x1, x2, max, eps1, eps2, &ind, snx);

	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;
}

/*********************************************************/
/* secant法（はさみうち法）による非線形方程式f(x)=0の解  */
/*      x1,x2 : 初期値（x1 < x2）                        */
/*      max : 最大試行回数                               */
/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
/*      eps2 : 終了条件２（｜f(x(k))｜＜eps2）           */
/*      ind : > 0 : 収束回数                             */
/*            =-1 : 収束しなかった                       */
/*      fun : f(x)を計算する関数名                       */
/*      return : 解                                      */
/*********************************************************/
#include <math.h>

double secant(double x1, double x2, int max, double eps1, double eps2,
              int *ind, double (*g)(double))
{
	double f, f1, f2, x = 0.0;
	int count;

	count = 0;
	*ind  = 0;
	f1    = (*g)(x1);
	f2    = (*g)(x2);

	if (fabs(f1) < eps2)
		x = x1;
	else {
		if (fabs(f2) < eps2)
			x = x2;
		else {
			while (*ind == 0) {
				count += 1;
				if (fabs(f2-f1) < eps2)
					*ind = -1;
				else {
					x = x2 - f2 * (x2 - x1) / (f2 - f1);
					f = (*g)(x);
					if (fabs(f) < eps2 ||
						 fabs(x2-x1) < eps1 || fabs(x2-x1) < eps1*fabs(x2))
						*ind = count;
					else {
						if (count < max) {
							if (f1*f2 < 0.0) {
								x2 = x;
								f2 = f;
							}
							else {
								x1 = x;
								f1 = f;
							}
						}
						else
							*ind = -1;
					}
				}
			}
		}
	}

	return x;
}