240.sagews

#
# RK4
#
# dy/dx = 0.2 * x * y
#

from pylab import *

# time step
h = 0.1

# initial condition
x0 = 1.0
y0 = 1.0

# final point
xn = 5.5

# 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
y = y0
print "xn   yn           e_abs        e_rel"
for x in arange(x0,xn+(h/2),h):
	k1 = fn(x, y)
	k2 = fn(x + 0.5 * h, y + 0.5 * h * k1)
	k3 = fn(x + 0.5 * h, y + 0.5 * h * k2)
	k4 = fn(x + h, y + h * k3)
	
        e_rel = abs(y - fn_exact(x))
	print "%.2lf" % x, "%e" % y, "%e" % e_rel, "%e" % ((e_rel / fn_exact(x)) * 100.0)
	y = y + (h/6) * (k1 + 2 * k2 + 2 * k3 + k4)

xn yn e_abs e_rel 1.00 1.000000e+00 0.000000e+00 0.000000e+00 1.10 1.021222e+00 5.252887e-11 5.143727e-09 1.20 1.044982e+00 1.312455e-10 1.255959e-08 1.30 1.071436e+00 2.471776e-10 2.306975e-08 1.40 1.100759e+00 4.151821e-10 3.771780e-08 1.50 1.133148e+00 6.550329e-10 5.780645e-08 1.60 1.168826e+00 9.927754e-10 8.493781e-08 1.70 1.208041e+00 1.462413e-09 1.210566e-07 1.80 1.251071e+00 2.108002e-09 1.684958e-07 1.90 1.298228e+00 2.986256e-09 2.300256e-07 2.00 1.349859e+00 4.169772e-09 3.089043e-07 2.10 1.406353e+00 5.751036e-09 4.089325e-07 2.20 1.468145e+00 7.847378e-09 5.345096e-07 2.30 1.535721e+00 1.060711e-08 6.906925e-07 2.40 1.609623e+00 1.421713e-08 8.832581e-07 2.50 1.690459e+00 1.891229e-08 1.118767e-06 2.60 1.778909e+00 2.498708e-08 1.404630e-06 2.70 1.875734e+00 3.280994e-08 1.749179e-06 2.80 1.981789e+00 4.284103e-08 2.161735e-06 2.90 2.098032e+00 5.565426e-08 2.652688e-06 3.00 2.225541e+00 7.196437e-08 3.233567e-06 3.10 2.365525e+00 9.266058e-08 3.917125e-06 3.20 2.519348e+00 1.188481e-07 4.717415e-06 3.30 2.688544e+00 1.518995e-07 5.649878e-06 3.40 2.874848e+00 1.935183e-07 6.731425e-06 3.50 3.080217e+00 2.458176e-07 7.980528e-06 3.60 3.306863e+00 3.114176e-07 9.417310e-06 3.70 3.557293e+00 3.935660e-07 1.106364e-05 3.80 3.834350e+00 4.962880e-07 1.294321e-05 3.90 4.141259e+00 6.245711e-07 1.508167e-05 4.00 4.481688e+00 7.845953e-07 1.750669e-05 4.10 4.859812e+00 9.840184e-07 2.024807e-05 4.20 5.280389e+00 1.232330e-06 2.333786e-05 4.30 5.748849e+00 1.541293e-06 2.681046e-05 4.40 6.271400e+00 1.925489e-06 3.070269e-05 4.50 6.855146e+00 2.403002e-06 3.505397e-05 4.60 7.508229e+00 2.996263e-06 3.990637e-05 4.70 8.239993e+00 3.733109e-06 4.530474e-05 4.80 9.061181e+00 4.648102e-06 5.129685e-05 4.90 9.984156e+00 5.784171e-06 5.793347e-05 5.00 1.102317e+01 7.194664e-06 6.526852e-05 5.10 1.219467e+01 8.945918e-06 7.335917e-05 5.20 1.351769e+01 1.112047e-05 8.226599e-05 5.30 1.501424e+01 1.382107e-05 9.205302e-05 5.40 1.670986e+01 1.717574e-05 1.027880e-04 5.50 1.863420e+01 2.134407e-05 1.145423e-04