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]

In [9]:

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

@author: LLP-admin
"""
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

def AdaptiveRK4(zk,vk,startstep,accuracy,toleranceFactor):
    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
avStep,avError=AdaptiveRK4(1,0,0.1,0.01,10)
print avStep
print avError

0.189250645105
0.000486186994821

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

Problem Set 3

Finite Difference for Ordinary Differential Equations

Problem 7.2

Problem 7.3

Problem 7.4