ルンゲ・クッタ法

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# ルンゲ・クッタ法  dx/dt=f(t,x)
#      time : 現在の時間
#      h : 時間刻み幅
#      x : 現在の状態
#      dx : 微係数（f(t,x):subで計算する）
#      g : 作業域（g[4][n]）
#      n : 微分方程式の次数
#      sub : 微係数を計算する関数の名前
#      return : time+h
#      coded by Y.Suganuma
############################################

def rkg(time, h, x, dx, g, n, sub) :

	h2 = 0.5 * h
	sub(time, x, dx)
	for i1 in range(0, n) :
		g[0][i1] = h * dx[i1]

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

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

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

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

	return time

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import rkg

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

			# 微係数の計算 */
def snx(time, x, dx) :
	dx[0] = x[1]
	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0
			# 初期設定
n    = 2
x    = np.array([0.0, 0.0], np.float)
dx   = np.empty(n, np.float)
g    = np.empty((4, n), np.float)
time = 0.0
h    = 0.01
			# 計算と出力
for i1 in range(0, 101) :
	print("time = " + str(time) + ", x = " + str(x[0]))
	time = rkg(time, h, x, dx, g, n, snx)