We may go a few steps further, before resorting to numerical means:
We start with Galileo's equation, and add a couple of initial conditions:
\ddot{\theta}+\omega^{2}\sin\theta=0,\theta(0)=\theta_{0},\dot{\theta}(0)=0
We multiply the equation with \dot{\theta} integrate, rearrange, and utilize intitial conditions, and gain:
\dot{\theta}=\pm\omega\sqrt{2(\cos\theta-\cos\theta_{0})}
The negative root is used on time intervals where \theta\to{-\theta_{0}}
(assuming \theta_{0}>0)
We thereby gain:
T=\frac{2}{\omega}\int_{-\theta_{0}}^{\theta_{0}}\frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_{0})}}
In the harmonic case, we have T_{h}=\frac{2\pi}{\omega}
Hence, given a relative error bound \epsilon we gain the bound of the initial angle as:
|1-\frac{1}{\pi}\int_{-\theta_{0}}^{\theta_{0}}\frac{d\theta}{\sqrt{2(\cos\theta-\cos\theta_{0})}}|<\epsilon