Reference request: error estimates for pendulum with small-angle approximation Let consider a simple pendulum of angular displacement $\theta(t)$ that starts at $\theta(0) = \alpha \in (-\pi,\pi)$ with no initial velocity. 
$\theta$ is given by the ODE
$$ \left\{ \begin{aligned}
& \theta^{\prime\prime}+\omega^2 \sin \theta = 0 \\
& \theta(0) = \alpha,\ \theta^\prime(0) = 0
\end{aligned} \right. $$
where $\omega := \sqrt{g/\ell}$, with $g$ the acceleration due to the gravity and $\ell$ is the length of the pendulum.
Physics textbooks tells us that under the small-angle approximation $\theta \ll 1$, since $\sin \theta \simeq \theta$, one gets the linearised equation
$$ \left\{ \begin{aligned}
& \theta_1^{\prime\prime}+\omega^2 \theta_1 = 0 \\
& \theta_1(0) = \alpha,\ \theta_1^\prime(0) = 0
\end{aligned} \right. $$
and one has $\theta \simeq \theta_1$.
Of course this is not very convincing... A lot of ODEs have chaotic behaviour, so it is not clear at all that replacing $\sin \theta$ by $\theta_1$ will not lead to huge changes to the solution.
I am looking for a rigorous mathematical analysis of this problem. I am particularly interested for a result that looks like an error estimate of the form
$$\forall t \geq 0,\quad \vert \theta(t) - \theta_1(t) \vert \leq F(t,\alpha)$$
where $F$ is as explicit as possible, with of course $F(t,\alpha) \to 0$ when $\alpha \to 0$.
I would also be interested for the high-order approximations $\theta_n$, where we replace $\sin \theta$ by the $n$ first terms of its Taylor series. For instance,
$$ \left\{ \begin{aligned}
& \theta_3^{\prime\prime}+\omega^2 \left( \theta_3 - \frac{1}{6} \theta_3^3+\frac{1}{150}\theta_3^5 \right) = 0 \\
& \theta_3(0) = \alpha,\ \theta_3^\prime(0) = 0
\end{aligned} \right. $$
Thanks!
 A: You won't get an error bound for all $t$ like you ask because the frequency of the approximation and the frequency of the real pendulum differ, so $\theta$ will drift away from the approximation.  The fact that the real solution is periodic comes from conservation of energy and the form of the equation.  The maximum amplitude has to repeat, then you are in the same state as the last maximum, so the swing will repeat, and the position is periodic.  
I don't believe you can get a closed form solution for the exact equation.  What we usually do is expand the solution in a series.  You write 
$$\theta''(t)+\omega^2 \theta(t)=\omega^2(\theta(t)-\sin(\theta(t))\approx \omega^2\left(\frac {\theta(t)^3}6-\frac{\theta(t)^5}{120}+\ldots\right)$$
The leading term cancels on the left, so you can solve the right ignoring the left, then plug each solution into the right and solve again.  You get a series expansion in harmonics of the original $\omega$
A: constant energy curves
As told in the comments, using the energy function simplifies the situation somewhat. Multiply with $2θ'(t)$ and integrate to get
$$
θ'(t)^2+2ω^2(1-\cosθ(t))=θ'(t)^2+\left(2ω\sin\frac{θ(t)}2\right)^2=2E=ω^2R^2=ω^2\left(2\sin\frac{α}2\right)^2
$$
It follows that the points $(\frac{θ'(t)}ω, 2\sin\frac{θ(t)}2)$ lie on a circle of the small radius $R$.
parametrization by modified angle
Now introduce $\phi$ as the angle on that circle, $\phi(0)=\frac\pi2$. Then
$$
θ'(t)=ωR\cosϕ(t), ~~2\sin\frac{θ(t)}2=R\sinϕ(t).
$$
Then also in the derivative of the second equation
$$
\cos\left(\frac{θ(t)}2\right)θ'(t)=R\cos(ϕ(t))ϕ'(t)
$$
connecting the pendulum angle and the modified phase space angle
Now insert the equation for $θ(t)$ and cancel common factors. This would only be problematic if $\cosϕ(t)=0$ constantly, that is,$θ'(t)=0$, a constant solution which could only happen for $α=0$.
After cancellation also eliminate $θ(t)$, then separate the variables
$$
ω\cos\left(\frac{θ(t)}2\right)=ϕ'(t)
\\
\implies ω=\left(1-\frac{R^2}{4}\sin^2ϕ(t)\right)^{-\frac12}ϕ'(t)
=ϕ'(t)+\frac{R^2}{8}\sin^2(ϕ(t))ϕ'(t)-\frac{R^4}{32}\sin^4(ϕ(t))ϕ'(t)\pm...
$$
Using $\sin^2ϕ=\frac{1-\cos(2ϕ)}2$, $\sin^4ϕ=\frac{\cos(4ϕ)-4\cos(2ϕ)+6}8$ etc., the integration of the right side will give trend terms and oscillation terms. With the two terms here this results in
\begin{align}
ωt&=\left[ϕ(t)-\frac\pi2\right]+\frac{R^2}{16}\left[ϕ(t)-\frac\pi2\right]-\frac{R^2}{32}\sin(2ϕ(t))-\frac{3R^4}{128}\left[ϕ(t)-\frac\pi2\right]+\frac{R^4}{128}\sin(2ϕ(t))-\frac{R^4}{1024}\sin(4ϕ(t))
\\
&=\left(1+\frac{R^2}{16}-\frac{3R^4}{128}+O(R^6)\right)\left[ϕ(t)-\frac\pi2\right]-\left(\frac{R^2}{32}-\frac{R^4}{128}+O(R^6)\right)\sin(2ϕ(t))
\\&~~~~ -\left(\frac{R^4}{1024}+O(R^6)\right)\sin(4ϕ(t))+O(R^6)
\end{align}
deriving approximation formulas
The linear approximation corresponds to setting $R=0$ in the last formula, so that $\phi(t)=\frac\pi2+ωt$. In the next order approximation one gets
$$
ϕ(t)=\frac\pi2+\left(1-\frac{R^2}{16}\right)ωt-\frac{R^2}{32}\sin(2ωt)
$$
so that one gets a correction to the frequency and at the same time small oscillations in the phase. The first implies a continuously growing phase difference in the small angle approximation, leading to a long period oscillation of the error of period $T=\frac{32\pi}{ωR^2}$,  containing about $16R^{-2}$ oscillations of the pendulum.
For the rather large value of $θ(0)=α=0.4$ I get then the following plot of numerical solution of the pendulum, the small angles approximation and the second order approximation, and the error of the approximations against the numerical solution.

As can be seen, the second order approximation gives a very good impression of what the error of the small angle approximation looks like, for a range of about $t< 40R^{-4}$
Using the second order approximation to replace the more exact numerical solution gives in the next plot a small but increasing phase shift in the first error plot in the next picture

A: Just a note on this subject.   
Starting from
$$
\theta '' + \omega ^{\,2} \sin \theta  = 0\quad \left| {\; - \pi  < \theta  \le \pi } \right.
$$
it is interesting to see what happens if we take half of the angle
$$
0 = 2\alpha '' + \omega ^{\,2} \sin 2\alpha  = 2\left( {\alpha '' + \left( {\omega ^{\,2} \cos \alpha } \right)\sin \alpha } \right)\quad \left| {\; - \pi /2 < \alpha  \le \pi /2} \right.
$$
So we have halved the oscillation, in change of introducing a variable frequency ($\cos \alpha$ is non negative).
Continuing we have
$$
0 = n\alpha '' + \omega ^{\,2} \sin n\alpha  = n\left( {\alpha '' + \left( {{1 \over n}\omega ^{\,2} \,U_{\,n - 1} (\cos \alpha )} \right)\sin \alpha } \right)\quad \left| {\; - \pi /n < \alpha  \le \pi /n} \right.
$$
where $U_{\,n - 1}$ is the Chebyshev polynomial of 2nd kind.
