# 7.2 - Harmonic Oscillator Approximations¶

In [64]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook

A=np.array([[0, 1],[-1, 0]])

def eulerUpdate(y, h, x):
return y + h*A.dot(y)

def rungeKuttaUpdate(y, h, x):
k1 = h*A.dot(y)
k2 = h*A.dot(y+k1/2.)
k3 = h*A.dot(y+k2/2.)
k4 = h*A.dot(y+k3)
return y+k1/6.+k2/3.+k3/3.+k4/6.

def adaptiveUpdate(y, h, x, maxError, r = 2., N=3):
error = float('inf')
difThreshold = 0.1*maxError
i=1

while True:
f1 = rungeKuttaUpdate(y, h, x)
f2_1 = rungeKuttaUpdate(y, h/2, x)
f2 = rungeKuttaUpdate(f2_1, h/2, x)
error = (f1-f2)[0]
if(abs(error)>maxError):
f3 = rungeKuttaUpdate(rungeKuttaUpdate(y, h/4, x),h/4, x)
error2 = (f2_1-f3)[0]
else:
if (maxError-abs(error))<difThreshold:
return (f1, h, error)

f2 = rungeKuttaUpdate(f1, h, x)
f3 = rungeKuttaUpdate(y, 2*h, x)
error2 = (f3-f2)[0]

if i>=N:
if abs(error2)>maxError:
return (f1, h, error)
else:
return (f2, 2*h, error)

h = h + min(max((abs(error)-maxError)*-gradient*r/i**0.5, -0.95*h),20*h)
i+=1

def localError(y, h, x, f):
f_1 = f(y,h,x)
f_2 = f(f(y, h/2, x), h/2,x)
val = f_1 - f_2
return val[0]

def error(y, x):
return y - np.array([np.cos(x), -np.sin(x)])

def calc(h, update, length=100*np.pi, initial = np.array([1,0])):
vals = np.ndarray((4, int(np.ceil(length/h))))
vals[0, :] = np.arange(0, length, h)
vals[1:3, 0] = initial
vals[3] = 0
x = 0
for i in range(1,vals.shape[1]):
y = vals[1:3, i-1]
vals[1:3,i] = update(y, h, x)
vals[3, i] = localError(y, h, x, update)
x+=h
return vals

def euler(h):
return calc(h, eulerUpdate)

def rungeKutta(h):
return calc(h, rungeKuttaUpdate)

h = initialH
vals = np.ndarray((5, 1))
vals[0, 0] = 0
vals[1:3, 0] = np.array([1, 0])
vals[3] = 0
x = 0
i = 1
while x < length:
y = vals[1:3, i-1]
val = np.ndarray(5)
y, h, error = adaptiveUpdate(y, h, x, maxError)
x+=h
val[0] = x
val[1:3] = y
val[3] = error
val[4] = h
vals = np.append(vals, val.reshape((5,1)), axis=1)
i+=1
return vals

def calcPlot(h):
valsEuler = euler(h)
valsRunge = rungeKutta(h)
plt.figure(figsize=(12, 8))
plt.title('Harmonic Oscillation Approximations with h=%f' % h)
cosRange = np.arange(0, 100*np.pi, 0.01)
cosVal = np.cos(cosRange)
plt.plot(cosRange, cosVal, label='System')
plt.plot(valsRunge[0,:], valsRunge[1,:], label='Runge-Kutta approximation')
plt.plot(valsEuler[0,:], valsEuler[1,:], label='Euler approximation')
plt.legend()
plt.ylim((-2,2))
plt.xlim((0,100*np.pi))
plt.show()

eulerError = error(valsEuler[1:3,-1], valsEuler[0,-1])
rungeError = error(valsRunge[1:3,-1], valsRunge[0,-1])

print(("Euler Error\n"
"=====\n"
"Mean abs local error:\t%f\n"
"Final value error:\t%f\n"
"Final slope error:\t%f\n"
"\n"
"Runge-Kutta Error\n"
"=====\n"
"Mean abs local error:\t%f\n"
"Final value error:\t%f\n"
"Final slope error:\t%f"
) % (np.mean(np.absolute(valsEuler[3,:])), eulerError[0], eulerError[1], np.mean(np.absolute(valsRunge[3,:])),rungeError[0],rungeError[1])
)

#plt.figure(figsize=(12, 8))
#plt.title('Errors over Harmonic Oscillation with Adaptive Approximations')
#plt.plot(valsRunge[0,:],valsRunge[3,:], label='local Error')
#plt.legend()
#plt.show()

plt.figure(figsize=(12, 8))
plt.title('Harmonic Oscillation Adaptive Approximations with maxError=%f' % maxError)
cosRange = np.arange(0, 100*np.pi, 0.01)
cosVal = np.cos(cosRange)
plt.plot(cosRange, cosVal, label='System')
plt.legend()
plt.ylim((-2,2))
plt.xlim((0,100*np.pi))
plt.show()

rungeError = error(vals[1:3,-1], vals[0,-1])

print(("Runge-Kutta Error\n"
"=====\n"
"Mean abs local error:\t%f\n"
"Final value error:\t%f\n"
"Final slope error:\t%f\n"
"Mean h:\t%f\n"
) % (np.mean(np.absolute(vals[3,:])),rungeError[0],rungeError[1], np.mean(vals[4,:]))
)

plt.figure(figsize=(12, 8))
plt.title('h over Harmonic Oscillation with Adaptive Approximations')
plt.plot(vals[0,:],vals[4,:], label='h')
plt.ylim((0,1))
plt.xlim((0,100*np.pi))
plt.legend()
plt.show()

plt.figure(figsize=(12, 8))
plt.title('Local Error Magnitudes over Harmonic Oscillation with Adaptive Approximations')
plt.plot((0, 100*np.pi), (maxError, maxError), 'red', label='maxError')
plt.plot(vals[0,:],np.absolute(vals[3,:]), label='local Error')
plt.ylim((0,2*maxError))
plt.xlim((0,100*np.pi))
plt.legend()
plt.show()

In [20]:
calcPlot(0.0001)

Euler Error
=====
Mean abs local error:	0.000000
Final value error:	0.015832
Final slope error:	0.000002

Runge-Kutta Error
=====
Mean abs local error:	0.000000
Final value error:	0.000000
Final slope error:	-0.000000

In [19]:
calcPlot(0.01)

Euler Error
=====
Mean abs local error:	0.000039
Final value error:	3.808983
Final slope error:	0.085658

Runge-Kutta Error
=====
Mean abs local error:	0.000000
Final value error:	-0.000000
Final slope error:	0.000000

In [18]:
calcPlot(1)

Euler Error
=====
Mean abs local error:	773281289593070220783632114572445969538875392.000000
Final value error:	-0.987344
Final slope error:	-182687704666362864775460604089535377456991567872.000000

Runge-Kutta Error
=====
Mean abs local error:	0.269047
Final value error:	-1.036285
Final slope error:	-0.020345


# 7.3 - Adaptive Harmonic Oscillation Approximations¶

In [21]:
calcAdaptivePlot(0.001)