ODE
7.4) Numerically solve the differential equation found in Problem 5.3:
l θ'' + (g + z'') sin(θ) = 0
z = A sin(w t)
Take the Motion of the platform to be periodic, and interactively explore the dynamics of the pendulum as a fuction of the amplitude and frequency of the excitation.
In [0]:
import math
from matplotlib import pyplot as plt
import numpy as np
In [0]:
## parameters
l = 10
g = 9.8
A = 5
w = 1.7
# initial time is 0
x = 0
# initial y(x), starting from 45 degree angle
y = 0.1
# initial dy(x) = 0, which means the initially the pendulum is stoped
dy = 0
In [0]:
## let θ to be y , t to be x
def current_ddy(x,y):
ddy = -1*(g - A*(w*w)*math.sin(w*x))*math.sin(y)/l
return ddy
In [0]:
def next_dy(x,y,h):
next_dy = dy + h*current_ddy(x,y)
return next_dy
In [0]:
def next_y(x,y,h):
next_y = y + h * dy
return next_y
In [0]:
## Location of Pendulum
def pLocation(x,y,length):
location = []
location.append(length*math.sin(y))
location.append(-1*length*math.cos(y))
location.append(A*math.sin(w*x))
return location
In [0]:
history_y = []
history_dy = []
history_ddy = []
history_x = []
history_px = []
ddy = current_ddy(x,y)
history_y.append(y)
history_dy.append(dy)
history_ddy.append(ddy)
history_x.append(x)
history_px.append(pLocation(x,y,l)[0])
In [8]:
h = 0.01
while ( x < 50000):
y = next_y(x,y,h)
dy = next_dy(x,y,h)
x = x + h
ddy = current_ddy(x,y)
history_y.append(y)
history_dy.append(dy)
history_dy.append(ddy)
history_x.append(x)
pLocation(x,y,l)
history_px.append(pLocation(x,y,l)[0])
print(y)
In [9]:
plt.plot(history_x,history_y,'b',c='red',label="theta")
#plt.plot(history_x,history_ddy,'b',c='red',label="theta")
plt.xlabel('time')
plt.ylabel('theta')
plt.legend()
plt.show()
