Period of a simple pendulum

Question

I'm writing a short program to simulate a simple pendulum. The equations of motion are

$$\frac{d\theta}{dt}=\omega$$ $$\frac{d\omega}{dt}=-\frac{g}{r}\sin\theta$$

For some small time step $dt$ the time evolution of $\theta$ and $\omega$ are governed by

$$\theta\left(t+dt\right)=\theta\left(t\right)+\omega\left(t\right)dt$$ $$\omega\left(t+dt\right)=\omega\left(t\right)-\frac{g}{r}\sin\left[\theta\left(t+dt\right)\right]dt$$

Assuming a small initial amplitude, the period of oscillation is

$$T=2\pi\sqrt{\frac{r}{g}}$$

Here is the python code I came up with, which works well

from numpy import *

file = open('task15.txt', 'w')
file.write("time" + "\t" + "theta" + "\t" + "omega" + "\n")

t = 0.0
tmax = 20.0
dt = 0.1
theta = 0.1
omega = 0.0
while (t <= tmax):

    theta += omega * dt
    omega -= sin(theta) * dt     #choose units such that (g/r) = 1
    file.write(str(t) + "\t" + str(theta) + "\t" + str(omega) "\n")
    t += dt

file.close()
exit()

By plotting the amplitude as a function of time, I have found graphically that the period asymptotically increases from the expected value (given above) as the initial amplitude increases, because the small angle approximation no longer holds.

I would like to update the code so it can calculate the period, although I'm not sure how to do this explicitly.

Somehow I need to create a user-defined function which returns the time when $\theta$ changes sign. Can I do this using the while loop? Then the difference between two successive times is half the period.

Can someone push me in the right direction with an algorithm and/or code? Thanks in advance.

Answer 1

If you want to include the period, I will recommend you compute it separately as it is much faster and accurate than your simulation.

Expanding the comment by semiclassical, the period of an arbitary amplitude simple pendulum has the form

$$T = 4 \sqrt{\frac{\ell}{2g}}\int_0^{\theta_0} \frac{d\theta}{\sqrt{\cos\theta - \cos\theta_0}} = 4\sqrt{\frac{\ell}{2g}} K\left(\sin\frac{\theta_0}{2}\right) $$ where $$K(k)\;\stackrel{def}{=}\; \int_0^{\pi/2} \frac{1}{\sqrt{1-k^2\sin^2 u}} du$$

is the complete elliptic integral of the first kind. However, there isn't any real problem to compute the integral on RHS numerically. This is because

$$K(k) = \frac{\pi}{2\text{AGM}(1-k,1+k)}\tag{*1}$$

where $\text{AGM}(x,y)$ is the Arithmetic-geometric mean of $x$ and $y$. To compute $\text{AGM}(x,y)$, you construct a pair of sequences $(x_n)$, $(y_n)$ by following rules:

$$(x_0,y_0) = (x,y)\quad\text{ and }\quad \begin{cases} x_n &= \frac12( x_{n-1}+y_{n-1} )\\ y_n &= \sqrt{x_{n-1}y_{n-1}} \end{cases} \quad\text{ for } n \ge 1 $$

The two sequences $(x_n)$, $(y_n)$ will converges to a common limit, the Arithmetic-geometric mean $\text{AGM}(x,y)$ we seek. The convergence is very fast, the "correct" number of digits grows exponentially with each iteration!

To compute $K(k)$ and hence the period, you just need to start from a pair of numbers $1-k$, $1+k$. Repeating taking AM and GM a few times to desired accuracy and plug it into $(*1)$.

Answer 2

1. Define a new variable "Tperiod" that you initialize to the value 0.

2. Because you detect the end of the period by the sign of "omega" you have to declare another variable omega2 that you put in at the end of your Loop. You can write omega2=omega to save the last value of Omega in another value.

3. In your loop "while (t <= tmax)" you insert the following if-condition:

"if (omega2/abs(omega2)<0 AND omega/abs(omega) > 0) then Tperiod = t".

Here the function $\frac{\omega}{|\omega|}=Omega/abs(omega)$ is the Signum function.

Answer 3

Only for small amplitudes (simple pendulum) you are allowed to take

$$T=2\pi\sqrt{\frac{r}{g}}$$

Time period T elapses between two consecutive times strictly when $ \theta=0 $ for this non-linear pendulum. A routine for second order ODE needs to be implemented.