Period of a simple pendulum 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.

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)$.
