%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()
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,"%"
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
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,"%"
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]
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]
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]