１変数関数の極小値求解プログラム２

/*		brent.c		*/
#include <stdio.h>
#include <math.h>

#define		ITMAX	100			/* 反復回数の上限 */
#define		CGOLD	0.3819660	/* 黄金分割比 */
#define		ZEPS	1.0e-10		/* 極小がちょうど x=0 にあるときは相対精度 tol */
								/* の代わりにこれを絶対精度とする */
#define		SHFT(a,b,c,d)	(a)=(b); (b)=(c); (c)=(d);
#define		SIGN(a, b)		((b) >= 0.0? fabs(a): -fabs(a))

/*
関数 f と、極小を囲い込む３点の座標 ax,bx,cx とを与えると、brent の方法で極小
を求める。与える３点の条件は、bx が ax,cx の間にあり、f(bx) が f(ax),f(bx) の
どちらよりも小さいことである。戻り値は極小値で、そのｘ座標は xmin に入る。
xmin の相対精度はほぼ tol になる。
*/
double brent(double ax, double bx, double cx, double (*f)(double),
		double tol, double *xmin)
{
	int iter;
	double a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm;
	double e = 0.0;						/* 前々回の更新量 */

	if(ax < cx)							/* a < b にする */
	{
		a = ax;
		b = cx;
	}
	else
	{
		a = cx;
		b = ax;
	}
	x = w = v = bx;						/* 初期化 */
	fw = fv = fx = (*f)(x);
	for(iter = 1; iter <= ITMAX; iter++)	/* 主ループ */
	{
		xm = 0.5 * (a + b);
		tol1 = tol * fabs(x) + ZEPS;
		tol2 = 2.0 * tol1;
		if(fabs(x - xm) <= (tol2 - 0.5 * (b - a)))	/* 収束判定 */
		{
			*xmin = x;					/* 最良の値を返す */
			return fx;
		}
		if(fabs(e) > tol1)				/* 放物線補間してみる */
		{
			r = (x - w) * (fx - fv);
			q = (x - v) * (fx - fw);
			p = (x - v) * q - (x - w) * r;
			q = 2.0 * (q - r);
			if(q > 0.0)	p = -p;
			q = fabs(q);
			etemp = e;
			e = d;
			if(fabs(p) >= fabs(0.5 * q * etemp)
				|| p <= q * (a - x)
				|| p >= q * (b - x))	/* 放物線補間の適否の検査 */
			{
				e = (x >= xm)? a - x: b - x;
				d = CGOLD * e;			/* 放物線補間は不適。大きい方の区間を */
										/* 黄金分割 */
			}
			else
			{
				d = p / q;				/* 放物線補間を採択する */
				u = x + d;
				if(u - a < tol2 || b - u < tol2)	d = SIGN(tol1, xm - x);
			}
		}
		else
		{
			e = (x >= xm)? a - x: b - x;
			d = CGOLD * e;
		}
		u = x + ((fabs(d) >= tol1)? d: SIGN(tol1, d));
		fu = (*f)(u);					/* 主ループでの関数値評価はここだけ */
		if(fu <= fx)
		{
			if(u >= x)	a = x;
			else		b = x;
			SHFT(v, w, x, u);
			SHFT(fv, fw, fx, fu);
		}
		else
		{
			if(u < x)	a = u;
			else		b = u;
			if(fu <= fw || w == x)
			{
				v = w;
				w = u;
				fv = fw;
				fw = fu;
			}
			else if(fu <= fv || v == x || v == w)
			{
				v = u;
				fv = fu;
			}
		}
	}
	fprintf(stderr, "Error : Too many iterations in brent.\n");
	*xmin = x;
	return fx;
}

double f(double x)
{
	return 5. + fabs((x - 1.6)*(x - 10.0)*(x + 8.0));
}

int main(void)
{
	double ax, bx, cx, fmin, tol, xmin;

	ax = 1.0;
	bx = 1.5;
	cx = 2.0;
	tol = 1.e-8;
	printf("f(x)=|(x-1.6)(x-10)(x+8)|\n");
	printf("極小点 x=1.6,10,-8\n");
	printf("[ input ]\n");
	printf("start point : ax=%f, bx=%f, cx=%f\n", ax, bx, cx);
	printf("tolerance = %e\n", tol);

	fmin = brent(ax, bx, cx, f, tol, &xmin);

	printf("[ output ]\n");
	printf("xmin=%f, fmin=%f\n", xmin, fmin);
	return 1;
}

f(x)=|(x-1.6)(x-10)(x+8)|
極小点 x=1.6,10,-8
[ input ]
start point : ax=1.000000, bx=1.500000, cx=2.000000
tolerance = 1.000000e-08
[ output ]
xmin=1.600000, fmin=0.000000

f(x)=|(x-1.6)(x-10)(x+8)|
極小点 x=1.6,10,-8
[ input ]
start point : ax=-10.000000, bx=0.000000, cx=20.000000
tolerance = 1.000000e-08
[ output ]
xmin=1.600000, fmin=5.000000

f(x)=|(x-1.6)(x-10)(x+8)|
極小点 x=1.6,10,-8
[ input ]
start point : ax=-10.000000, bx=-7.000000, cx=0.000000
tolerance = 1.000000e-08
[ output ]
xmin=-8.000000, fmin=5.000007

f(x)=|(x-1.6)(x-10)(x+8)|
極小点 x=1.6,10,-8
[ input ]
start point : ax=5.000000, bx=8.000000, cx=12.000000
tolerance = 1.000000e-08
[ output ]
xmin=10.000000, fmin=5.000007