オイラー法

今回はオイラー法。例として解いているのは dy/dx = (x^2+x+1)-(2x+1)y+y^2。
某教科書と同じものを使った。
次のサイトが参考になる。

・3 オイラー法
・オイラー法

def euler(x0, xn, y0, n):
	x = x0
	y = y0
	h = (xn - x0) / n    #刻み幅を求める。今回は0.1になる。
	ans = [[x, y]]
	for j in range(n):
		y += f(x, y) * h
		x += h
		tmp = [x, y]
		ans.append(tmp)
	return ans

def f(x, y):    #解く微分方程式
	return ((x**2 + x + 1) - (2*x + 1) * y + y**2)

def output(l, f):
	for i, j in l:
		f.write('%f %f\n' % (i, j))

def main():
	n = 20    #分割数
	x0 = 0.0    #xの初期値
	xn = 2.0    #終了地点
	y0 = 0.5    #yの初期値
	ans = euler(x0, xn, y0, n)
	for i in ans:
		print i
	f = open('euler.dat', 'w')
	output(ans, f)
	f.close()

if __name__ == '__main__':
	main()

結果はこちら

>>> 
[0.0, 0.5]
[0.1, 0.575]
[0.2, 0.6500625]
[0.30000000000000004, 0.725311875390625]
[0.4, 0.8008697069863916]
[0.5, 0.876852388485688]
[0.6, 0.953368921907856]
[0.7, 1.0305189892141025]
[0.7999999999999999, 1.1083913705158035]
[0.8999999999999999, 1.1870627572050847]
[0.9999999999999999, 1.2665969841419948]
[1.0999999999999999, 1.347044680923156]
[1.2, 1.4284433202680829]
[1.3, 1.510817623298785]
[1.4000000000000001, 1.5941802679982415]
[1.5000000000000002, 1.6785328388464043]
[1.6000000000000003, 1.7638669524164194]
[1.7000000000000004, 1.850165494984202]
[1.8000000000000005, 1.937403913074167]
[1.9000000000000006, 2.0255515052995596]
[2.0000000000000004, 2.114572672817902]

さらにgnuplotで描いてみるとこんな感じになる。

元の関数というか厳密解のグラフはこんな感じ。

描いている範囲が違うからこの2枚をそのまま重ねても意味ないけど、実際に比較してみると結構誤差がある。