239.sagews

#
# Comparison of Euler's method and RK4
#

from pylab import *

# time step
h = 0.1

# initial condition
x0 = 1.0
y0 = 1.0

# final point
xn = 15

# function (DE)
def fn(x, y):
	return 0.2 * x * y;

# exact solution
def fn_exact(x):
	return exp(0.1*(x*x - 1))

# algorithm
yrk = y0
yeu = y0

xr=arange(x0,xn,h)
yr = []			# ys for rk4
ye = []			# ys for euler
yx = []			# ys for exact
for x in xr:
	# RK4
	k1 = fn(x, yrk)
	k2 = fn(x + 0.5 * h, yrk + 0.5 * h * k1)
	k3 = fn(x + 0.5 * h, yrk + 0.5 * h * k2)
	k4 = fn(x + h, yrk + h * k3)
	yr.append(yrk)
	yrk = yrk + (h/6) * (k1 + 2 * k2 + 2 * k3 + k4)
	# Euler's
	ye.append(yeu)
	yeu = yeu + h * fn(x, yeu)
	# exact
	yx.append(fn_exact(x))

# show me!
figure()
plot(xr,yr, 'b', label="rk4")	# rk4 is blue
plot(xr,ye, 'r', label="euler")	# euler is red
plot(xr,yx, 'g', label="exact")	# exact is green
legend()

savefig("rk-vs-euler.png")