非線形方程式(ニュートン法)

<?php

/****************************/
/* ニュートン法             */
/* (exp(x)-3.0*x=0の根)   */
/*      coded by Y.Suganuma */
/****************************/

	$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));

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

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

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

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

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

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

	return $x;
}

?>