# Problem Set 3

Problem 7.1:

## Problem 7.2

In [8]:
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 06 04:39:21 2014

@author: Marwan Sarieddine
"""

from math import *

#Exact solution y=cosx
def Euler(step):
t=0
vk=0 #y'=0
zk=1 #y =1
zkold=zk
vold=vk
yEuler=[]
errorLocal=[]
averageLocal=0.0
errorFinal=0.0
SlopeFinal=0.0
ytrue=[]
while t<100*pi:
yEuler.append(zkold)
t+=step
vnew=vold-zkold*step
zknew=zkold+vold*step
vold=vnew
zkold=zknew
errorLocal.append(abs(zkold-cos(t)))
averageLocal=sum(errorLocal)/int((100*pi/step))
errorFinal=abs(zkold-cos(100*pi))
SlopeFinal=vold
return t,yEuler,ytrue,averageLocal,errorFinal,SlopeFinal,step

Table=[["Step","Average Local Error","Final Error   ", "Final Slope Error"]]
for i in range(5):
x,yEuler,ytrue,averageLocal,errorFinal,SlopeFinal,step=Euler(1.0/10**i)
Table.append([step,averageLocal,errorFinal,SlopeFinal])
for i in range(len(Table)):
print Table[i]

x=[1,3]
print sum(x)

def RK4(step):
zkold=1.0
vold=0.0
errorLocal=[]
averageLocal=0.0
errorFinal=0.0
SlopeFinal=0.0
t=0
while t < 100*pi:
t+=step
v1=-1*step*zkold
z1=step*vold
v2=-1*step*(zkold+z1/2.0)
z2=step*(vold+v1/2.0)
v3=-1*step*(zkold+z2/2.0)
z3=step*(vold+v2/2.0)
v4=-1*step*(zkold+z3)
z4=step*(vold+v3)
vnew=vold+v1/6.0+v2/3.0+v3/3.0+v4/6.0
zknew=zkold+z1/6.0+z2/3.0+z3/3.0+z4/6.0
vold=vnew
zkold=zknew
errorLocal.append(abs(zkold-cos(t)))

averageLocal=sum(errorLocal)/int((100*pi/step))
errorFinal=abs(zkold-cos(100*pi))
SlopeFinal=vold
return averageLocal,errorFinal,SlopeFinal,step

Table=[["Step Size","   Average Local Error   ","     Final Error   ", "Final Slope"]]

for i in range(5):
averageLocal,errorFinal,SlopeFinal,step=RK4(1.0/10**i)
Table.append([step,averageLocal,errorFinal,SlopeFinal])
for i in range(len(Table)):
print Table[i]

['Step', 'Average Local Error', 'Final Error   ', 'Final Slope Error']
[1.0, 1.551487937718581e+45, 1.8268770466636286e+47, -1.8268770466636286e+47]
[0.1, 247040.06672304886, 3321022.606325886, 5176276.8086821195]
[0.01, 0.9076946141525466, 3.809889267410002, 0.04683407654122684]
[0.001, 0.05272470123543376, 0.17008889324935872, -0.00073706376090179]
[0.0001, 0.00502628740211624, 0.015831982938209643, -3.412567986018993e-05]
4
['Step Size', '   Average Local Error   ', '     Final Error   ', 'Final Slope']
[1.0, 0.4583836958232742, 0.9113033309273068, 0.11566814182709298]
[0.1, 8.338424087197775e-05, 0.000840724185820596, -0.040461811176643964]
[0.01, 8.314401574435761e-09, 2.700476563477494e-07, -0.000734614775384406]
[0.001, 2.0247207335593923e-10, 2.6984871737134597e-07, -0.0007346409519864463]
[0.0001, 2.091999950367963e-09, 5.99935323641887e-10, -3.464102085443581e-05]



## Problem 7.3

In [9]:
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 10 17:13:42 2014

"""
from math import *

def RK4(zk,vk,step):
zkold=zk
vold=vk
t=0
while t < 100*pi:
t+=step
v1=-1*step*zkold
z1=step*vold
v2=-1*step*(zkold+z1/2.0)
z2=step*(vold+v1/2.0)
v3=-1*step*(zkold+z2/2.0)
z3=step*(vold+v2/2.0)
v4=-1*step*(zkold+z3)
z4=step*(vold+v3)
vnew=vold+v1/6.0+v2/3.0+v3/3.0+v4/6.0
zknew=zkold+z1/6.0+z2/3.0+z3/3.0+z4/6.0
vold=vnew
zkold=zknew
return zkold,vold

zkold=zk
vold=vk
t=0
step=startstep
error=0
stepTotal=0
count=0
while t<100*pi:
y1,v1=RK4(zkold,vold,step)
yhalf,vhalf=RK4(zkold,vold,step/2.0)
ysecondhalf,vsecondhalf=RK4(yhalf,vhalf,step/2.0)
Delta=ysecondhalf-y1
if Delta>accuracy:
step=step*(abs(accuracy/Delta))**0.2
elif Delta<accuracy/toleranceFactor:
step=step*(abs(accuracy/Delta))**0.2
else:
zkold=y1
vold=v1
error=abs(cos(t)-zkold)
stepTotal+=step
t+=step
count+=1
averageStep=stepTotal/count
averageError=error/count
return  averageStep,averageError
print avStep
print avError

0.189250645105
0.000486186994821



## Problem 7.4

In [13]:
import matplotlib.pylab as plt
import numpy as np
from math import *
%pylab inline
count=0
step=0.001
tolerance=0.0001
vk=1 #y'=0
zk=0 #y =1
zkold=zk
vold=vk
thetaEuler=[]
#y''+y=0 <=> v'=-y
# v=y'
#vk+1=vk-y(step)
#zk+1=zk+vk(step)
error=100
g=9.81
l=20.0
A=2.0
w=4
for i in range(20000):
thetaEuler.append(zkold)
count+=1
vnew=vold-(sin(zkold)*(g/l+A*sin(w*i*step)))*step
zknew=zkold+vold*step
vold=vnew
zkold=zknew
total= count*step

t=np.linspace(0,total,count)

plt.plot(t,thetaEuler,'r')
plt.show()

Populating the interactive namespace from numpy and matplotlib


WARNING: pylab import has clobbered these variables: ['plt', 'sinh', 'trunc', 'step', 'tan', 'gamma', 'cosh', 'radians', 'sin', 'fmod', 'expm1', 'ldexp', 'exp', 'frexp', 'ceil', 'copysign', 'isnan', 'cos', 'degrees', 'tanh', 'fabs', 'sqrt', 'floor', 'hypot', 'log', 'log10', 'e', 'pi', 'log1p', 'modf', 'isinf']
%pylab --no-import-all prevents importing * from pylab and numpy


In []: