7.2 Simple Harmonic Oscillator

$$ \ddot{y} + y = 0 $$

initial conditions:

$$ y(0) = 1 $$$$ \dot{y(0)} = 0 $$

Find $y(t)$ from $t=0$ to $100\pi$

analytical solution:

$$y = cos(t)$$

import numpy as np
import math
import matplotlib.pyplot as plt

dt = .2

timeSteps = np.arange(0, 100*math.pi, dt)
positions = np.zeros(len(timeSteps))
velocities = np.zeros(len(timeSteps))
error = np.zeros(len(timeSteps))

#initial conditions
positions[0] = 1
velocities[0] = 0

def exactSol(t):
    return math.cos(t)

def f(currentY, currentV):
    return [currentV, -currentY]

def stepEuler(currentY, currentV, _dt):
    nextAccel = -currentY
    nextVelocity = nextAccel*_dt + currentV
    nextPosition = nextVelocity*_dt + currentY
    return [nextPosition, nextVelocity]

def stepRK4(currentY, currentV, _dt):
    k1 = _dt*np.array(f(currentY, currentV))
    k2 = _dt*np.array(f(currentY + k1[0]/2, currentV + k1[1]/2))
    k3 = _dt*np.array(f(currentY + k2[0]/2, currentV + k2[1]/2))
    k4 = _dt*np.array(f(currentY + k3[0], currentV + k3[1]))
    return (np.array([currentY, currentV]) + k1/6 + k2/3 + k3/3 + k4/6).tolist()

for i in range(len(positions)):
    if (i == 0): continue
    nextVals = stepEuler(positions[i-1], velocities[i-1], dt)
    positions[i] = nextVals[0]
    velocities[i] = nextVals[1]
    error[i] = abs(exactSol(timeSteps[i]) - positions[i])
    
plt.plot(timeSteps, positions)
plt.plot(timeSteps, error, 'r')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Euler (local error in red)')
axes = plt.gca()
axes.set_xlim([0,314])
axes.set_ylim([-1.5,1.5])
plt.show()

for i in range(len(positions)):
    if (i == 0): continue
    nextVals = stepRK4(positions[i-1], velocities[i-1], dt)
    positions[i] = nextVals[0]
    velocities[i] = nextVals[1]
    error[i] = abs(exactSol(timeSteps[i]) - positions[i])
    
plt.plot(timeSteps, positions)
plt.plot(timeSteps, error, 'r')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('RK4 - constant step')
axes = plt.gca()
axes.set_xlim([0,314])
axes.set_ylim([-1.5,1.5])
plt.show()

/usr/local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

def calcEuler(_dt):
    timeSteps = np.arange(0, 100*math.pi, _dt)
    positions = np.zeros(len(timeSteps))
    velocities = np.zeros(len(timeSteps))
    error = np.zeros(len(timeSteps))
    
    #initial conditions
    positions[0] = 1
    velocities[0] = 0
    for i in range(len(positions)):
        if (i == 0): continue
        nextVals = stepEuler(positions[i-1], velocities[i-1], _dt)
        positions[i] = nextVals[0]
        velocities[i] = nextVals[1]
        error[i] = abs(exactSol(timeSteps[i]) - positions[i])
    return [np.mean(np.array(error)), abs(exactSol(timeSteps[i]) - positions[i])]#avg local error, final error

numSteps = 10
localError = np.zeros(numSteps)
finalError = np.zeros(numSteps)
stepSizes = np.zeros(numSteps)
index = 0
for i in np.linspace(0.001, 1, numSteps):
    errors = calcEuler(i)
    stepSizes[index] = i
    localError[index] = errors[0]
    finalError[index] = errors[1]
    index = index+1
    
plt.plot(stepSizes, localError, label='local error')
plt.plot(stepSizes, finalError, label='final error')
plt.ylabel('error')
plt.xlabel('step size')
plt.title('Euler error')
axes = plt.gca()
axes.set_xlim([0,1])
axes.set_ylim([0,2])
plt.legend()
plt.show()

def calcRK4(_dt):
    timeSteps = np.arange(0, 100*math.pi, _dt)
    positions = np.zeros(len(timeSteps))
    velocities = np.zeros(len(timeSteps))
    error = np.zeros(len(timeSteps))
    
    #initial conditions
    positions[0] = 1
    velocities[0] = 0
    for i in range(len(positions)):
        if (i == 0): continue
        nextVals = stepRK4(positions[i-1], velocities[i-1], _dt)
        positions[i] = nextVals[0]
        velocities[i] = nextVals[1]
        error[i] = abs(exactSol(timeSteps[i]) - positions[i])
    return [np.mean(np.array(error)), abs(exactSol(timeSteps[i]) - positions[i])]#avg local error, final error

numSteps = 10
localError = np.zeros(numSteps)
finalError = np.zeros(numSteps)
stepSizes = np.zeros(numSteps)
index = 0
for i in np.linspace(0.001, 1, numSteps):
    errors = calcRK4(i)
    stepSizes[index] = i
    localError[index] = errors[0]
    finalError[index] = errors[1]
    index = index+1
    
plt.plot(stepSizes, localError, label='local error')
plt.plot(stepSizes, finalError, label='final error')
plt.ylabel('error')
plt.xlabel('step size')
plt.title('RK4 error')
axes = plt.gca()
axes.set_xlim([0,1])
axes.set_ylim([0,2])
plt.legend()
plt.show()

7.3 fourth order RK4 adaptive stepper

time = []
positions = []
velocities = []
error = []
stepSizes = []

def exactSol(t):
    return math.cos(t)

def calcRK4Adaptive(_dt, _desError):
    
    #initial conditions
    positions.append(1)
    velocities.append(0)
    time.append(0)
    error.append(0);
    stepSizes.append(0);
    
    t = _dt
    
    while t<math.pi*100:
        
        lookingForStepSol = True
        i = len(positions)
        
        exact = exactSol(t)
                
        while lookingForStepSol:
            
            #full step
            fullStep = stepRK4(positions[i-1], velocities[i-1], _dt)

            #two half steps
            halfStep = stepRK4(positions[i-1], velocities[i-1], _dt/2)
            halfStep = stepRK4(halfStep[0], halfStep[1], _dt/2)

            errorEst = abs(fullStep[0] - halfStep[0])
        
#             print _dt
            if (errorEst > _desError):
                 _dt *= 1/1.3#decrease the step size and try again
            else :
                stepSizes.append(_dt)
                _dt *= 1.2#increase the step size
                lookingForStepSol = False
           
        positions.append(fullStep[0])
        velocities.append(fullStep[1])
        error.append(abs(exactSol(t) - positions[i]))
        time.append(t)
        t += _dt

desError = 0.0000000001
calcRK4Adaptive(0.001, desError)
# print error

plt.plot(time, positions, label='rk4 solution')
linTime = np.linspace(0, math.pi*100, 1000)
plt.plot(linTime, np.cos(linTime), label='exact')
plt.ylabel('error')
plt.xlabel('step size')
plt.title('Adaptive RK4')
axes = plt.gca()
axes.set_xlim([0,100*math.pi])
axes.set_ylim([-1.5,1.5])
plt.legend()
plt.show()

Numerically solve:

$$l\ddot{\theta} + ( g + \ddot{z} )sin\theta = 0$$

motion of the platform is periodic, interactively explore the dynamics of the pendulum as a func of amp and freq of excitation

g = 9.8
l = 1

_dt = 0.1

timeSteps = np.arange(0, 100, _dt)
positions = np.zeros(len(timeSteps))
velocities = np.zeros(len(timeSteps))
error = np.zeros(len(timeSteps))

def plot_pend(platAmp, platFreq):

    timeSteps = np.arange(0, 100, _dt)
    positions = np.zeros(len(timeSteps))
    velocities = np.zeros(len(timeSteps))
    error = np.zeros(len(timeSteps))

    #initial conditions
    positions[0] = 0.1
    velocities[0] = 0
    for i in range(len(positions)):
        if (i == 0): continue
        #z = platAmp*sin(2*pi*platFreq*t)
        zAccel = platAmp*(2*math.pi*platFreq)**2*math.sin(2*math.pi*platFreq*timeSteps[i])
        nextAccel = -(g+zAccel)/l*math.sin(positions[i-1])
        nextVelocity = nextAccel*_dt + velocities[i-1]
        nextPosition = nextVelocity*_dt + positions[i-1]
        velocities[i] = nextVelocity
        positions[i] = nextPosition
        error[i] = abs(exactSol(timeSteps[i]) - positions[i])

    plt.plot(timeSteps, positions, label='theta: amp %.2f, freq: %.2f' %(platAmp, platFreq))
    # plt.plot(timeSteps, platAmp*np.sin(2*math.pi*platFreq*np.array(timeSteps)), label='platform (z)')
    # plt.ylabel('theta')
    plt.xlabel('time')
    plt.title('Pend on Platform')
    axes = plt.gca()
    axes.set_xlim([0,100])
    # axes.set_ylim([0,2])
    

plot_pend(0.1, 5)
axes = plt.gca()
axes.set_xlim([0,100])
plt.legend()
plt.show()

plot_pend(1, 1)
axes = plt.gca()
axes.set_xlim([0,100])
plt.legend()
plt.show()

plot_pend(10, 2)
axes = plt.gca()
axes.set_xlim([0,100])
plt.legend()
plt.show()