Using Euler's method to compute the frequency of a nonlinear pendulum In my studies of numerical methods I have come across the following exercise:

We consider the following second-order ODE $$\ddot{\theta}+\sin(\theta) = 0 $$ and we reduce it to a two-dimensional system of first-order ODEs using the 2-vector $$y(t) = \begin{pmatrix} y_1(t) \\ y_2(t) \end{pmatrix} = \begin{pmatrix} \theta(t) \\ \dot{\theta}(t) \end{pmatrix} $$ to get the first-order 2-dimensional system of ODEs $$ 
y'(t) = \begin{pmatrix} y_2(t) \\ -\sin(y_1(t)) \end{pmatrix} = f(y(t)) $$ We are asked to use Euler's method given by the recurrence relation $$ y_n = y_{n-1} + hf(y_{n-1}) $$
A.We are asked to choose initial conditions for $ y(0) $and use Euler's method to compute one full swing of the pendulum.


B. We are asked to use part A with a step size $h$ small enough to compute the frequency of the pendulum to get 3 digits of accuracy and explain.


C. We are asked to compute the theoretical error bound (as hints, we are asked to look at the Lipschitz constant of $f(y)$, to use energy conservation to obtain a bound on $\lVert y''(t) \rVert $ and to turn the error estimate for $y$ into the error estimate for the period of one swing).

I did part A with a computer, but I do not know how to do parts B and C. For one thing, part B has me stumped. When my stepsize $h$ is very small on the order of $10^{-3}$ I see a periodic solution which I can extract a time period from $T$ and the frequency would be $ f=\frac{1}{T} $. But here is my question: how can I find and justify a stepsize $h$ small enough so that $f$ is computed with 3 digits of accuracy? This is what has me stumped. I think I also need help on part C. I thank all helpers.
 A: B. In first order, you get the root $T_h$ of a function $f_h(t)=a(t)+hb(t)$ as approximation of a root $T_*$ of $a$. Expanding $f_h(T_*+hv)=a(T_*)+\dot a(T_*)hv+hb(T_*)+O(h^2)$ for $T_h=T_*+hv$ gives an estimate of $v=-\frac{ b(T_*)}{\dot a(T_*)}$.

the numerical approximation $f_h(t)=θ_h(t)$ of step size $h$ seen as $a(t)+hb(t)+...$ with $a(t)=θ_0(t)$ the exact solution, for several values of $h$. It is visible that not only the vertical value has a perturbation proportional to $h$, but also the root location.
So if you know an approximation for $\dot a(T_h)$ and $b(T_h)$, you get $T_h+h\frac{b(T_h)}{\dot a(T_h)}$ as improved root estimation of the root $T_*$ of $a=θ_0$. The important part is that $h\frac{b(T_h)}{\dot a(T_h)}$ is an error estimate of $T_h$.
$\dot a(T_h)=\dot θ_0(T_h)$ you get directly from the differential equation, $b(T_h)$ can be estimated by comparing the results for two different step sizes.
Which raises the question of if it is simpler to just estimating the error of $T_h$ by comparing it to $T_{2h}$. So compute the error for some relatively large but still reasonable $h$ and then scale $h$ so that the expected scaled-down error is in the desired region.

Numerical values for $T_h$ with a secant line. The slope for small $h$ is a little smaller than $0.5$, but still this rough estimate is sufficient to determine $h=10^{-3}$ as sufficient to get 3 correct digits after the dot.
\begin{array}{c|c}
h&T_h\\\hline
0.000500  & 6.70013638\\
0.001000 &  6.70029805\\
0.002000  & 6.70062424
\end{array}
C. just asks for bounds on $|b(T)|$ based on the global error formula $$e(T)\le\frac{M_2}{2L}(e^{LT}-1)h$$ where $M_2$ is a bound for the second derivative around the solution and $L$ the Lipschitz constant. $\dot a(T)$ can again be used directly, if you are looking for a root of $a(t)=\dot θ(t)$, then the value of $\dot a(t)=\ddot θ(t)=-\sinθ(t)$ is known approximately because $θ(T)$ will still be close to the maximum amplitude.
