ルンゲ・クッタ法（4次）

修正オイラー法と同じ要領でテイラー展開した式を4次の項まで残してやると精度があがる。詳しくは次のサイトで。

・ルンゲクッタ法
・ルンゲ＝クッタ法 - Wikipedia
・微分方程式の数値計算（ルンゲ・クッタ法）-数学アルゴリズム演習ノート-

def runge_kutta(x0, y0, dx, endx):
	x = x0
	y = y0
	res = [[x, y]]
	while x <= endx:
		s1 = f(x, y)
		s2 = f(x + dx/2, y + dx*s1/2)
		s3 = f(x + dx/2, y + dx*s2/2)
		s4 = f(x + dx, y + dx*s3)
		x += dx
		y += (s1 + 2*s2 + 2*s3 + s4) * dx / 6
		res.append([x, y])
	return res

def f(x, y):
	return (x**2 +x + 1) - (2*x + 1)*y + y**2

def output(filename, ans):
	f = open(filename, 'w')
	for i in ans:
		f.write('%f %f\n' % (i[0], i[1]))
	f.close()

def main():
	x0 = 0.0    #初期値
	y0 = 0.5
	dx = 0.1    #刻み幅
	endx = 2.0    #終了地点
	filename = 'runge_kutta.txt'
	ans = runge_kutta(x0, y0, dx, endx)
	output(filename, ans)

if __name__ == '__main__':
	main()

ファイルの中身はこちら。

0.000000 0.500000
0.100000 0.575021
0.200000 0.650166
0.300000 0.725557
0.400000 0.801312
0.500000 0.877541
0.600000 0.954344
0.700000 1.031812
0.800000 1.110026
0.900000 1.189051
1.000000 1.268941
1.100000 1.349740
1.200000 1.431475
1.300000 1.514165
1.400000 1.597816
1.500000 1.682426
1.600000 1.767982
1.700000 1.854465
1.800000 1.941851
1.900000 2.030109
2.000000 2.119203

他の方法との比較は後で気が向いたらやることにする。