非線形方程式(セカント法)

/*************************************/
/* セカント法による 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.secant(x1, x2, max, eps1, eps2, ind);

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

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

abstract class Secant
{
	/*********************************************************/
	/* secant法(はさみうち法)による非線形方程式f(x)=0の解  */
	/*      x1,x2 : 初期値(x1 < x2)                        */
	/*      max : 最大試行回数                               */
	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
	/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
	/*      ind : > 0 : 収束回数                             */
	/*            =-1 : 収束しなかった                       */
	/*      return : 解                                      */
	/*********************************************************/

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

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

		count  = 0;
		ind[0] = 0;
		f1     = snx(x1);
		f2     = snx(x2);

		if (Math.abs(f1) < eps2)
			x = x1;
		else {
			if (Math.abs(f2) < eps2)
				x = x2;
			else {
				while (ind[0] == 0) {
					count += 1;
					if (Math.abs(f2-f1) < eps2)
						ind[0] = -1;
					else {
						x = x2 - f2 * (x2 - x1) / (f2 - f1);
						f = snx(x);
						if (Math.abs(f) < eps2 ||
						    Math.abs(x2-x1) < eps1 || Math.abs(x2-x1) < eps1*Math.abs(x2))
							ind[0] = count;
						else {
							if (count < max) {
								if (f1*f2 < 0.0) {
									x2 = x;
									f2 = f;
								}
								else {
									x1 = x;
									f1 = f;
								}
							}
							else
								ind[0] = -1;
						}
					}
				}
			}
		}

		return x;
	}
}