微分方程式（ルンゲ・クッタ）

<?php

/*
          初期設定
*/
	$n  = 2;
	$x  = array($n);
	$dx = array($n);
	$g  = array(4);
	for ($i1 = 0; $i1 < 4; $i1++)
		$g[$i1] = array($n);
	$time = 0.0;
	$h    = 0.01;
	$x[0] = 0.0;
	$x[1] = 0.0;
/*
          計算と出力
*/
	for ($i1 = 0; $i1 < 101; $i1++) {
		printf("time = %f, x = %f\n", $time, $x[0]);
		$time = rkg($time, $h, $x, $dx, $g, $n, "snx");
	}

/****************/
/* 微係数の計算 */
/****************/
function snx($time, $x, &$dx)
{
	$dx[0] = $x[1];
	$dx[1] = -2.0 * $x[0] - 3.0 * $x[1] + 1.0;
}

/*******************************************/
/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
/*      time : 現在の時間                  */
/*      h : 時間刻み幅                     */
/*      x : 現在の状態                     */
/*      dx : 微係数（f(t,x):snxで計算する）*/
/*      g : 作業域（g[4][n]）              */
/*      n : 微分方程式の次数               */
/*      sub : 微係数を計算する関数の名前   */
/*      return : time+h                    */
/*******************************************/
function rkg($time, $h, &$x, $dx, $g, $n, $sub)
{
	$h2 = 0.5 * $h;

	$sub($time, $x, $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[0][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $x[$i1] + 0.5 * $g[0][$i1];
	$sub($time, $g[1], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $x[$i1] + 0.5 * $g[1][$i1];
	$sub($time, $g[2], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $x[$i1] + $g[2][$i1];
	$sub($time, $g[3], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$x[$i1] = $x[$i1] + ($g[0][$i1] + 2.0 * $g[1][$i1] + 2.0 * $g[2][$i1] + $g[3][$i1]) / 6.0;

	return $time;
}

?>