ルンゲ・クッタ法

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

#***************/
# 微係数の計算 */
#***************/
snx = Proc.new { |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 : 微分方程式の次数               */
#      snx : 微係数を計算する関数の名前   */
#      return : time+h                    */
#******************************************/
def rkg(time, h, x, dx, g, n, &sub)

	h2 = 0.5 * h

	sub.call(time, x, dx)
	for i1 in 0 ... n
		g[0][i1] = h * dx[i1]
	end

	time += h2
	for i1 in 0 ... n
		g[1][i1] = x[i1] + 0.5 * g[0][i1]
	end
	sub.call(time, g[1], dx)
	for i1 in 0 ... n
		g[1][i1] = h * dx[i1]
	end

	for i1 in 0 ... n
		g[2][i1] = x[i1] + 0.5 * g[1][i1]
	end
	sub.call(time, g[2], dx)
	for i1 in 0 ... n
		g[2][i1] = h * dx[i1]
	end

	time += h2
	for i1 in 0 ... n
		g[3][i1] = x[i1] + g[2][i1]
	end
	sub.call(time, g[3], dx)
	for i1 in 0 ... n
		g[3][i1] = h * dx[i1]
	end

	for i1 in 0 ... n
		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0
	end

	return time
end

		# 初期設定
n  = 2
x  = Array.new(n)
dx = Array.new(n)
g  = Array.new(4)
for i1 in 0 ... 4
	g[i1] = Array.new(n)
end
time = 0.0
h    = 0.01
x[0] = 0.0
x[1] = 0.0
		# 計算と出力
for i1 in 0 ... 101
	printf("time = %f, x = %f\n", time, x[0])
	time = rkg(time, h, x, dx, g, n, &snx)
end