非線形方程式(二分法)

/*********************************/
/* 二分法による exp(x)-3x=0 の根 */
/*      coded by Y.Suganuma      */
/*********************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double eps1, eps2, x, x1, x2;
		int max, ind[] = new int [1];
	/*
	          データの設定
	*/
		eps1 = 1.0e-10;
		eps2 = 1.0e-10;
		max  = 100;
		x1   = 0.0;
		x2   = 1.0;
	/*
	          実行と結果
	*/
		Kansu kn = new Kansu();
		x = kn.bisection(x1, x2, eps1, eps2, max, ind);

		System.out.println("   ind=" + ind[0] + "  x=" + x + "  f=" + kn.snx(x));
	}
}

/****************/
/* 関数値の計算 */
/****************/
class Kansu extends Bisection
{
	double snx(double x)
	{
		double y = Math.exp(x) - 3.0 * x;
		return y;
	}
}

abstract class Bisection
{
	/*********************************************************/
	/* 二分法による非線形方程式(f(x)=0)の解                  */
	/*      x1,x2 : 初期値                                   */
	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
	/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
	/*      max : 最大試行回数                               */
	/*      ind : 実際の試行回数                             */
	/*            (負の時は解を得ることができなかった)     */
	/*      return : 解                                      */
	/*********************************************************/

	abstract double snx(double x);   // 定義しておく必要あり

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

		f1 = snx(x1);
		f2 = snx(x2);

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

		else {
			ind[0] = 0;
			if (f1*f2 == 0.0)
				x0 = (f1 == 0.0) ? x1 : x2;
			else {
				sw = 0;
				while (sw == 0 && ind[0] >= 0) {
					sw      = 1;
					ind[0] += 1;
					x0      = 0.5 * (x1 + x2);
					f0      = snx(x0);
					if (Math.abs(f0) > eps2) {
						if (ind[0] <= max) {
							if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
								sw = 0;
								if (f0*f1 < 0.0) {
									x2 = x0;
									f2 = f0;
								}
								else {
									x1 = x0;
									f1 = f0;
								}
							}
						}
						else
							ind[0] = -1;
					}
				}
			}
		}

		return x0;
	}
}