Chapter-26 : Stiffness & Multi step Methods

Ex:26.1 Pg: 754

In [4]:
%matplotlib inline
from matplotlib.pyplot import plot,xlabel,ylabel,show,title,legend
from numpy import arange,exp
def f(t,y):
    yp=-1000*y+3000-2000*exp(-t)
    return yp
y0=0#
#explicit euler
h1=0.0005#
y1 =[y0]
count=1#
t=arange(0,0.0061,0.0001)
for i in arange(0,0.00591,0.0001):
    y1.append(y1[(count-1)]+f(i,y1[(count-1)])*h1)
    count=count+1#

plot(t,y1)
h2=0.0015#
y2=[y0]#
count=1#
t=arange(0,0.0061,0.0001)
for i in arange(0,0.00591,0.0001):
    y2.append(y2[(count-1)]+f(i,y2[(count-1)])*h2)
    count=count+1#

plot(t,y2)
title("y vs t")
xlabel("t")
ylabel("y")
h=legend(["h-0.0005","h=0.0015"])
show()
#implicit order
h3=0.05#
y3=[y0]#
count=1#
t=arange(0,0.401,0.01)
for j in arange(0,0.40,0.01):
    y3.append((y3[(count-1)]+3000*h3-2000*h3*exp(-(j+0.01)))/(1+1000*h3))
    count=count+1#

plot(t,y3)
title("y vs t")
xlabel("t")
ylabel("y")
show()

Ex:26.2 Pg: 758

In [5]:
from math import exp
print "f(x,y) = 4*exp(0.8*x) - 0.5*y"
#f'(x,y) = 4*exp(0.8*x) - 0.5*y
h = 1#
x=range(0,5)
y = [2]
x1 = -1#
y1 = -0.3929953#
y10 = y1 + (4*exp(0.8*x[0]) - 0.5*y[0])*2
y11 = y[0] + (4*exp(0.8*x[0]) - 0.5*y[0] + 4*exp(0.8*x[0]) - 0.5*y10)*h/2
y12 = y[0] + (3 + 4*exp(0.8*x[1]) - 0.5*y11)*h/2#
t = 6.360865#
y20 = y[0] + (4*exp(0.8*x[1]) - 0.5*t) *2
y21 = t + (4*exp(0.8*x[1]) - 0.5*t + 4*exp(0.8*x[2]) - 0.5*y20)*h/2
print "the first corrector yields y = ",y21
t = 14.84392
e = (t - y21)*100/t#
print "error = ",e,"%"
f(x,y) = 4*exp(0.8*x) - 0.5*y
the first corrector yields y =  15.7669298488
error =  -6.21810039929 %

Ex:26.3 Pg: 762

In [6]:
x1 = 1#
x2 = 2#
y1 = 6.194631#
y2 = 14.84392#
y10 = 5.607005#
y11 = 6.360865#
y20 = 13.44346#
y21 = 15.30224#
Ec1 = -(y11 - y10)/5#
print "Ec (x = 1) = ",Ec1
e1 = y1 - y11#
print "true error (x = 1) = ",e1
Ec2 = -(y21 - y20)/5#
print "Ec (x = 2) = ",Ec2
e2 = y2 - y21#
print "true error (x = 2) = ",e2
Ec (x = 1) =  -0.150772
true error (x = 1) =  -0.166234
Ec (x = 2) =  -0.371756
true error (x = 2) =  -0.45832

Ex:26.4 Pg: 763

In [7]:
x0 = 0#
x1 = 1#
x2 = 2#
y1 = 6.194631#
y2 = 14.84392#
y10 = 5.607005#
y11 = 6.360865#
y1m = y11 - (y11 - y10)/5#
e = (y1 - y1m)*100/y1#
print "ym =",y1m
print "error = ",e,"%"
y20 =2+(4*exp(0.8*x1) - 0.5*y1m)*2#
e2 = (y2 - y20)*100/y2#
print "y20 = ",y20
print "error = ",e2,"%"
y2o = y20 + 4* (y11 - y10)/5#
e2 = (y2 - y2o)*100/y2#
print "y20 = ",y2o
print "error = ",e2,"%"
y21 = 15.21178#
y23 = y21 - (y21 - y20)/5#
print "y2 = ",y23
e3 = (y2 - y23)*100/y2#
print "error = ",e3,"%"
ym = 6.210093
error =  -0.249603245133 %
y20 =  13.5942344279
error =  8.41883796235 %
y20 =  14.1973224279
error =  4.35597586123 %
y2 =  14.8882708856
error =  -0.2987814916 %

Ex:26.5 Pg: 773

In [9]:
from numpy import arange,exp
print "f(x,y) = 4*exp(0.8*x) - 0.5*y"
#f'(x,y) = 4*exp(0.8*x) - 0.5*y
h = 1#
x = arange(-3,4.1,h)
y=[-4.547302,-2.306160,-0.3929953,2,2]
y1=[0,0,0,0]
for  i in range(3,7):
    y.append(y[(i-3)] + 4*h*(2*(4*exp(0.8*x[(i)]) - 0.5*y[(i)]) - 4*exp(0.8*x[(i-1)]) + 0.5*y[(i-1)] + 2*(4*exp(0.8*x[(i-2)]) - 0.5*y[(i-2)]))/3)
    y1.append(y[(i-1)] + h*(4*exp(0.8*x[(i-1)]) - 0.5*y[(i-1)] +4 * (4*exp(0.8*x[(i)]) - 0.5*y[(i)]) + 4*exp(0.8*x[(i+1)]) - 0.5*y[(i+1)])/3)

print "x = ",x[4:8]
print "y0 = ",y[4:8]
print "corrected y1 = ",y1[4:8]
f(x,y) = 4*exp(0.8*x) - 0.5*y
x =  [ 1.  2.  3.  4.]
y0 =  [2, 6.0227228499844756, 20.083112587836688, 41.872835141608761]
corrected y1 =  [6.9056637734795867, 19.136474376155903, 38.385462010990274, 82.755745651588995]

Ex:26.6 Pg: 774

In [13]:
from numpy import arange,exp
print "f(x,y) = 4*exp(0.8*x) - 0.5*y"
#f'(x,y) = 4*exp(0.8*x) - 0.5*y
h = 1#
x = arange(-3,4.1,h)
y = [-4.547302,-2.306160,-0.3929953,2]
m= [0,0,0,0,y[3]]
for i in range(3,7):
    y.append(y[(i)] + h *(55* (4*exp(0.8*x[(i)]) - 0.5*y[(i)]) / 24 - 59 * (4*exp(0.8*x[(i-1)]) - 0.5*y[(i-1)]) / 24 + 37*(4*exp(0.8*x[(i-2)]) - 0.5*y[(i-2)])/24 - 9*(4*exp(0.8*x[(i-3)]) - 0.5*y[(i-3)])/24))
    m.append(y[(i+1)])
    y.append(y[(i)] + h*(9*(4*exp(0.8*x[(i+1)]) - 0.5*y[(i+1)])/24 +19*(4*exp(0.8*x[(i)]) - 0.5*y[(i)])/24 - 5*(4*exp(0.8*x[(i-1)]) - 0.5*y[(i-1)])/24 + (4*exp(0.8*x[(i-2)]) - 0.5*y[(i-2)])/24))

print "x = ",x[4:8]
print "y0 = ",m[4:8]
print "y1 = ",y[4:8]
f(x,y) = 4*exp(0.8*x) - 0.5*y
x =  [ 1.  2.  3.  4.]
y0 =  [2, 6.0075392692969114, 6.2532143855636217, 14.488238413222703]
y1 =  [6.0075392692969114, 6.2532143855636217, 14.488238413222703, 16.39224775082873]

Ex:26.7 Pg: 775

In [17]:
from numpy import arange,exp
#dy/dx = -y
#y = exp(-x)
h = 0.5#
x = arange(-1.5,10.1,h)
y=[exp(-x[0]),exp(-x[1]),exp(-x[2]),1]
m=[0,0,0,y[3]]
for  i in range(3,23):
    y.append(y[(i-3)] + 4*h*(2*(-y[(i)]) + y[(i-1)] + 2*(-y[(i-2)]))/3)
    m.append(y[(i+1)])
    y.append(y[(i-1)] + h*(-y[(i-1)] +4 * (-y[(i)]) + (-y[(i+1)]))/3)

print "\nx = ",x[3:23]
print "\ny0(milnes method) = ",m[3:23]
print "\ncorrected y1(milnes method) = ",y[3:23]
for i in range(3,23):
    y[(i+1)] = y[(i)] + h *(55* (-y[(i)]) / 24 - 59 * (-y[(i-1)]) / 24 + 37*(-y[(i-2)])/24 - 9*(-y[(i-3)])/24)#
    m[(i+1)] = y[(i+1)]
    y[(i+1)] = y[(i)] + h*(9*(-y[(i+1)])/24 +19*(-y[(i)])/24 - 5*(-y[(i-1)])/24 + (-y[(i-2)])/24)#

print "\ny0(fourth order adams method) = ",m[3:23]
print "\ny1(fourth order adams method) = ",y[3:23]
x =  [ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5  7.
  7.5  8.   8.5  9.   9.5]

y0(milnes method) =  [1, 0.62312747952608971, 0.60341314566242521, 0.35581682815742122, 0.3173461560388694, -0.073744603832379019, 0.05769466447055438, 0.097016353530268939, 0.21274204327392865, -0.36734037063735558, 0.09724068674400127, 0.43888095058828192, 0.3040024218422181, -0.1934006684100632, -0.1160863384290636, 0.32478026608871907, -0.052055689173705413, -0.36964931180484378, 0.00024232753216428538, 0.55995471612928904]

corrected y1(milnes method) =  [1, 0.62312747952608971, 0.60341314566242521, 0.35581682815742122, 0.3173461560388694, -0.073744603832379019, 0.05769466447055438, 0.097016353530268939, 0.21274204327392865, -0.36734037063735558, 0.09724068674400127, 0.43888095058828192, 0.3040024218422181, -0.1934006684100632, -0.1160863384290636, 0.32478026608871907, -0.052055689173705413, -0.36964931180484378, 0.00024232753216428538, 0.55995471612928904]

y0(fourth order adams method) =  [1, 0.27152768648678072, 0.65650943823969443, -0.097035945241965349, 0.90661897200241937, -0.41645854530303128, -0.006339328628468005, 0.0808211944939453, -0.031548281043396395, -0.01627399436000334, 0.0039276970657568626, -0.002208005756050492, -0.0027943948998587478, -0.00042783320198906672, -0.00032845203812321728, -0.00042142028381515442, -0.00017932875240256902, -8.560953717861569e-05, -7.1995388648036731e-05, -4.2020211202092437e-05]

y1(fourth order adams method) =  [1, 0.56419948638876194, 0.68342664355034821, -0.010131383289966767, -0.11667598191807985, -0.0076991693311919962, -0.015405634984349756, -0.02283279228603189, -0.0093238635936440419, -0.0046392587796136846, -0.0040348813308916029, -0.0023127486894963844, -0.0011969855369064737, -0.00079981133060367615, -0.00049813831748820084, -0.00028031874518611127, -0.00017096177282523818, -0.00010605960387014003, -6.2547086111287706e-05, -3.7396246727791117e-05]