非線形方程式（ニュートン法）

/****************************/
/* ニュートン法             */
/* （exp(x)-3.0*x=0の根）   */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

double newton(double(*)(double), double(*)(double), double, double,
              double, int, int *);
double snx(double);
double dsnx(double);

int main()
{
	double eps1, eps2, x, x0;
	int max, ind;

	eps1 = 1.0e-7;
	eps2 = 1.0e-10;
	max  = 20;
	x0   = 0.0;

	x = newton(snx, dsnx, x0, eps1, eps2, max, &ind);

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

	return 0;
}

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

double newton(double(*f)(double), double(*df)(double), double x0,
              double eps1, double eps2, int max, int *ind)
{
	double g, dg, x, x1;
	int sw;

	x1   = x0;
	x    = x1;
	*ind = 0;
	sw   = 0;

	while (sw == 0 && *ind >= 0) {

		sw    = 1;
		*ind += 1;
		g     = (*f)(x1);

		if (fabs(g) > eps2) {
			if (*ind <= max) {
				dg = (*df)(x1);
				if (fabs(dg) > eps2) {
					x = x1 - g / dg;
					if (fabs(x-x1) > eps1 && fabs(x-x1) > eps1*fabs(x)) {
						x1 = x;
						sw = 0;
					}
				}
				else
					*ind = -1;
			}
			else
				*ind = -1;
		}
	}

	return x;
}

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

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