非線形方程式(二分法)

<?php

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

/*
          データの設定
*/
	$eps1 = 1.0e-10;
	$eps2 = 1.0e-10;
	$max  = 100;
	$x1   = 0.0;
	$x2   = 1.0;
/*
          実行と結果
*/
	$x = bisection("snx", $x1, $x2, $eps1, $eps2, $max, $ind);

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

/****************/
/* 関数値の計算 */
/****************/
function snx($x)
{
	return exp($x) - 3.0 * $x;
}

/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解                  */
/*      fn : f(x)を計算する関数名                        */
/*      x1,x2 : 初期値                                   */
/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
/*      max : 最大試行回数                               */
/*      ind : 実際の試行回数                             */
/*            (負の時は解を得ることができなかった)     */
/*      return : 解                                      */
/*********************************************************/
function bisection($fn, $x1, $x2, $eps1, $eps2, $max, &$ind)
{
	$x0 = 0.0;
	$f1 = $fn($x1);
	$f2 = $fn($x2);

	if ($f1*$f2 > 0.0)
		$ind = -1;

	else {
		$ind = 0;
		if ($f1*$f2 == 0.0)
			$x0 = ($f1 == 0.0) ? $x1 : $x2;
		else {
			$sw = 0;
			while ($sw == 0 && $ind >= 0) {
				$sw    = 1;
				$ind += 1;
				$x0    = 0.5 * ($x1 + $x2);
				$f0    = $fn($x0);

				if (abs($f0) > $eps2) {
					if ($ind <= $max) {
						if (abs($x1-$x2) > $eps1 && abs($x1-$x2) > $eps1*abs($x2)) {
							$sw = 0;
							if ($f0*$f1 < 0.0) {
								$x2 = $x0;
								$f2 = $f0;
							}
							else {
								$x1 = $x0;
								$f1 = $f0;
							}
						}
					}
					else
						$ind = -1;
				}
			}
		}
	}

	return $x0;
}

?>