Not-so-simple harmonic motion

When the bob makes an angle $\theta$ with the downward vertical, conservation of energy tells us that $\tfrac12m\ell^2 \dot\theta^2 + mg\ell(1-\cos\theta) \; =\; mg\ell$ and hence $\dot\theta^2 \; = \; \frac{2g}{\ell}\cos\theta$ Alternatively, we could obtain this equation by considering the transverse equation of circular motion: $m\ell \ddot{\theta} \; = \; -mg\sin\theta$ and integrating it once. Thus the time for the bob to sweep from one side to the other is $\begin{array}{rcl} \displaystyle T & = & \displaystyle \sqrt{\frac{\ell}{2g}} \int_{-\frac12\pi}^{\frac12\pi} \frac{d\theta}{\sqrt{\cos\theta}} \; = \; \sqrt{\frac{2\ell}{g}} \int_0^{\frac12\pi} \frac{d\theta}{\sqrt{\cos\theta}} \; = \; \sqrt{\frac{2\ell}{g}} \int_0^{\frac12\pi} \frac{\cos\theta\,d\theta}{(1 - \sin^2\theta)^{\frac34}} \\ & = & \displaystyle \sqrt{\frac{2\ell}{g}} \int_0^1 \frac{dx}{(1-x^2)^{\frac34}} \; = \; \sqrt{\frac{\ell}{2g}} \int_0^1 \frac{dx}{x^{\frac12} (1-x)^{\frac34}} \; = \; \sqrt{\frac{\ell}{2g}}B\big(\tfrac12,\tfrac14\big) \end{array}$ which makes $T = \boxed{0.64227}$ .

This integral can also be evaluated in terms of elliptic integrals, but the beta function is probably more familiar!

Who needs to solve elliptic integrals when you've got your very own personal computer? I feel like I have cheated...

The equation of motion for the pendulum is:

$\displaystyle \ddot{\theta} = -\frac{g}{l} \sin(\theta)$

Here's a simulation I did using Runge-Kutta, and the results are strikingly accurate. The equation of motion is fed in, and is broken into two first order linked equations to apply the RK algorithm:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

import math

initialPosition = math.pi/2

time = 0
deltaT = 10**-7

length = 0.3
g = 10

u_1 = initialPosition
u_2 = 0

def differentialEquation(t,x,xdot):  #applying RK4 by splitting 2nd order ODE into 2 first order ODEs
    return -(g/length)*math.sin(x)

while u_1 >= -1.57079632679:
    #evaluating slope at various points

    k0 = deltaT*u_2
    L0 = deltaT*differentialEquation(time,u_1,u_2)

    k1 = deltaT*(u_2 + 0.5*L0)
    L1 = deltaT*differentialEquation(time + 0.5*deltaT, u_1 + 0.5*k0, u_2 + 0.5*L0)

    k2 = deltaT*(u_2 + 0.5*L1)
    L2 = deltaT*differentialEquation(time + 0.5*deltaT, u_1 + 0.5*k1, u_2 + 0.5*L1)

    k3 = deltaT*(u_2 + L2)
    L3 = deltaT*differentialEquation(time + deltaT, u_1 + k2, u_2 + L2)
    #RK4 approximation of values

    u_1 = u_1 + (1/6)*(k0 + 2*k1 + 2*k2 + k3)
    u_2 = u_2 + (1/6)*(L0 + 2*L1 + 2*L2 + L3)

    time += deltaT

print("time: " +str(time))

As the angle through which the pendulum has been released is not small enough we cannot treat its motion as SHM.So now the usual way of circular motion works.obtain an integral.this can be evaluated by application of maclaurin series to the denominator.the correction in the period is to be noted.this is a form of elliptic integrals.