Chapter-25 : Runge-Kutta Methods

Ex:25.1 Pg: 708

from numpy import arange
#dy/dx = -2*x**3 + 12*x**2 - 20*x + 8.5
#therefore, y = -0.5*x**4 + 4*x**3 - 10*x**2 + 8.5 + c
x1 = 0#
y1 = 1#
h = 0.5#
c =-(-0.5*x1**4 + 4*x1**3 - 10*x1**2 + 8.5*x1 - y1)#
x = arange(0,4.1,0.5)
print "x = ",x
y=[]
for xx in x:
    y.append(-0.5*xx**4 + 4*xx**3 - 10*xx**2 + 8.5*xx + c)
print "\ntrue values of y = ",y
fxy=[]
for xx in x:
    fxy.append(-2*xx**3 + 12*xx**2 - 20*xx + 8.5)
y2=[y[0]]
e = [(y[0] - y2[0]) * 100 / y[0]]
for i in range(1,9):
    y2.append(y2[(i-1)] + fxy[(i-1)]*h)
    e.append((y[(i)] - y2[(i)])*100/y[(i)])

print "\ny by euler method =",y2
print "\nerror =",e

x =  [ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4. ]

true values of y =  [1.0, 3.21875, 3.0, 2.21875, 2.0, 2.71875, 4.0, 4.71875, 3.0]

y by euler method = [1.0, 5.25, 5.875, 5.125, 4.5, 4.75, 5.875, 7.125, 7.0]

error = [0.0, -63.106796116504853, -95.833333333333329, -130.98591549295776, -125.0, -74.712643678160916, -46.875, -50.993377483443709, -133.33333333333334]

Ex:25.2 Pg: 712

#f(x,y) = dy/dx = -2*x**3 + 12*x**2 - 20*x + 8.5
#f'(x,y) = -6*x**2 + 24*x - 20
#f"(x,y) = -12*x + 24
#f"'(x,y) = -12
x = 0#
Et2 = (-6*x**2 + 24*x - 20) * 0.5**2 / 2
Et3 = (-12*x + 24) * (0.5)**3 / 6#
Et4 = (-12) *(0.5 ** 4) / 24#
Et = Et2 + Et3 + Et4#
print "Total truncation error =",Et

Total truncation error = -2.03125

Ex:25.3 Pg: 713

from numpy import arange
#dy/dx = -2*x**3 + 12*x**2 - 20*x + 8.5
#therefore, y = -0.5*x**4 + 4*x**3 - 10*x**2 + 8.5 + c
x1 = 0#
y1 = 1#
h = 0.25#
c =-(-0.5*x1**4 + 4*x1**3 - 10*x1**2 + 8.5*x1 - y1)#
xx = arange(0,4.1,h)
print "x = ",xx
y=[]
for x in xx:
    y.append(-0.5*x**4 + 4*x**3 - 10*x**2 + 8.5*x + c)
print "\ntrue values of y = ",y
fxy=[]
for x in xx:
    fxy.append(-2*x**3 + 12*x**2 - 20*x + 8.5)
y2= [y[0]]
e = [(y[0] - y2[0]) * 100 / y[0]]
for i in range(1,17):
    y2.append(y2[(i-1)] + fxy[(i-1)]*h)
    e.append((y[(i)] - y2[(i)])*100/y[(i)])

print "\ny by euler method =",y2
print "\nerror =",e

x =  [ 0.    0.25  0.5   0.75  1.    1.25  1.5   1.75  2.    2.25  2.5   2.75
  3.    3.25  3.5   3.75  4.  ]

true values of y =  [1.0, 2.560546875, 3.21875, 3.279296875, 3.0, 2.591796875, 2.21875, 1.998046875, 2.0, 2.248046875, 2.71875, 3.341796875, 4.0, 4.529296875, 4.71875, 4.310546875, 3.0]

y by euler method = [1.0, 3.125, 4.1796875, 4.4921875, 4.34375, 3.96875, 3.5546875, 3.2421875, 3.125, 3.25, 3.6171875, 4.1796875, 4.84375, 5.46875, 5.8671875, 5.8046875, 5.0]

error = [0.0, -22.04424103737605, -29.854368932038835, -36.986301369863014, -44.791666666666664, -53.127354935945739, -60.2112676056338, -62.267839687194524, -56.25, -44.569939183318851, -33.045977011494251, -25.073056691992985, -21.09375, -20.741699008193187, -24.337748344370862, -34.662437698232893, -66.666666666666671]