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

/**********************************/
/* 微分方程式（ルンゲ・クッタ法） */
/*      例：d2x/dt+3dx/dt+2x=1    */
/*      coded by Y.Suganuma       */
/**********************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double time, h, x[], dx[], g[][];
		int i1, n = 2;
	/*
	          初期設定
	*/
		x    = new double [n];
		dx   = new double [n];
		g    = new double [4][n];

		time = 0.0;
		h    = 0.01;
		x[0] = 0.0;
		x[1] = 0.0;
	/*
	          計算と出力
	*/
		Kansu kn = new Kansu();

		for (i1 = 0; i1 < 101; i1++) {
			System.out.println("time = " + time + " x = " + x[0]);
			time = kn.rkg(time, h, x, dx, g, n);
		}
	}
}

/****************/
/* 微係数の計算 */
/****************/
class Kansu extends RKG
{
	void snx(double time, double x[], double dx[])
	{
		dx[0] = x[1];
		dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
	}
}

abstract class RKG
{
	/*************************************************/
	/* ルンゲ・クッタ法  dx/dt=f(t,x)                */
	/*      time : 現在の時間                        */
	/*      h : 時間刻み幅                           */
	/*      x : 現在の状態                           */
	/*      dx : 微係数（f(t,x):snxで計算する）      */
	/*      g : 作業域（g[4][n]）                    */
	/*      n : 微分方程式の次数                     */
	/*      return : time+h                          */
	/*************************************************/

	abstract void snx(double time, double x[], double dx[]);

	double rkg(double time, double h, double x[], double dx[], double g[][], int n)
	{
		int i1;
		double h2;

		h2 = 0.5 * h;

		snx(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];
		snx(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];
		snx(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];
		snx(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;
	}
}