1. Linear solution#

The linear solution is the first order Taylor series approximation for the equation of motion

$\ddot{\theta }(t)=-\frac{g}{L}\sin \theta \approx -\frac{g}{L}\theta$

The solution to this simple harmonic oscillator function is

$\theta (t)={\theta }_0\cos \omega t+\frac{{\dot{\theta }}_0}{\omega }\sin \omega t$

where ${\theta }_0$ is the initial angle, ${\dot{\theta }}_0$ is the initial angular velocity [in rad/s], and $\omega =\sqrt{\frac{g}{L}}$.

Here, you can plot the solution for

• $L=1m$

• $g=9.81m/s^2$

• ${\theta }_0=\frac{\pi }{3}=[{60}^o]$

• ${\dot{\theta }}_0=0rad/s$

g = 9.81 # m/s/s
L = 1 # m
w = np.sqrt(g/L) # rad/s
t = np.linspace(0, 4*np.pi/w) # 2 time periods of motion
theta0 = np.pi/3
dtheta0 = 0

theta = theta0*np.cos(w*t) + dtheta0/w*np.sin(w*t)

plt.plot(t, theta*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')

Text(0, 0.5, '$\\theta(t)$ [degrees]')

In the plot above, you should see harmonic motion. If it starts from rest at ${60}^o$ the pendulum will swing to $-{60}^o$ in $\frac{T}{2}=\frac{\pi }{\omega }$ seconds.

2. Numerical solution#

In the numerical solution, there is no need to linearize the equation of motion. You do need to create a state variable as such

$\mathbf{x}=[x_1,x_2]=[\theta ,\dot{\theta }]$

now, you can describe 2 first order equations in a function as such

1. ${\dot{x}}_1=x_2$

2. ${\dot{x}}_2=\ddot{\theta }=-\frac{g}{L}\sin \theta$

putting this in a Python function it becomes

g = 9.81
L = 1
def pendulum(t, x):
    '''pendulum equations of motion for theta and dtheta/dt
    arguments
    ---------
    t: current time
    x: current state variable [theta, dtheta/dt]
    outputs
    -------
    dx: current derivative of state variable [dtheta/dt, ddtheta/ddt]'''
    
    dx = np.zeros(len(x))
    dx[0] = x[1]
    dx[1] = -g/L*np.sin(x[0])
    return dx

This function pendulum defines the differential equation. Next, integrate with solve_ivp to find $\theta (t)$ as such.

Note: Here, I am using some previously defined values:

• theta0: the initial angle

• dtheta0: the initial angular velocity

• t: the timesteps from the previous linear analysis $t=(0...2T)$

Show code cell source Hide code cell source

from scipy.integrate import solve_ivp

sol = solve_ivp(pendulum, [0, t.max()], [theta0, 0], t_eval = t)
plt.plot(sol.t, sol.y[0]*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')

Text(0, 0.5, '$\\theta(t)$ [degrees]')

Comparing results#

At first glance, it would seem both solutions are close enough, but here you can plot them together and compare.

plt.plot(t, theta*180/np.pi)
plt.plot(sol.t, sol.y[0]*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')

Text(0, 0.5, '$\\theta(t)$ [degrees]')

You should notice a difference in time period between the two solutions, which means the natural frequencies are different.

You can also consider angular velocity and angular acceleration, below a set of 3 plots are shown for $\theta (t),\dot{\theta }(t),and\ddot{\theta }(t)$.

dtheta = -theta0*w*np.sin(w*t) + dtheta0*np.cos(w*t)
ddtheta = -theta0*w**2*np.cos(w*t) - dtheta0*w*np.sin(w*t)

ddtheta_numerical = -g/L*np.sin(sol.y[0]) # using equation of motion

plt.figure(figsize= (10, 11))
plt.subplot(311)
plt.plot(t, theta*180/np.pi)
plt.plot(sol.t, sol.y[0]*180/np.pi)

plt.ylabel(r'$\theta(t)$ [degrees]')

plt.subplot(312)
plt.plot(sol.t, sol.y[1])
plt.plot(t, dtheta)

plt.ylabel(r'$\dot{\theta}(t)$ (rad/s)')


plt.subplot(313)
plt.plot(t, ddtheta)
plt.plot(t, ddtheta_numerical)
plt.xlabel('time (s)')
plt.ylabel(r'$\ddot{\theta}(t)$ (rad/s$^2$)')

Text(0, 0.5, '$\\ddot{\\theta}(t)$ (rad/s$^2$)')