Simple but accurate periodic solutions for the nonlinear pendulum equation

Despite its elementary structure, the simple pendulum oscillations are described by a nonlinear differential equation whose exact solution for the angular displacement from vertical as a function of time cannot be expressed in terms of an elementary function, so either a numerical treatment or some analytical approximation is ultimately demanded. Such solutions have been thoroughly investigated due to the abundance of distinct pendular systems in nature and, more recently, due to the availability of automatic data acquisition systems in undergraduate laboratories. However, it is well-known that numerical solutions to differential equations usually loose accuracy (due to accumulation of roundoff errors) and polynomial approximations diverge after long time intervals. In this work, I take a few terms of the Fourier series expansion of the elliptic function sn(u; k) as a source of accurate periodic solutions for the pendulum equation. Interestingly, these approximations remain accurate for arbitrarily long time intervals, even for large amplitudes, which shows its adequacy for the analysis of experimental data gathered in classical mechanics classes.


Introduction
The simple pendulum is a prototype for the study of nonlinear systems and their stability [1].In fact, pendular oscillations arise in distinct fields of physical science and technology, e.g.acoustic vibrations [2], molecular oscillations [3], optically torqued nanorods [4], Josephson superconducting junctions [3,5], elliptic filters in electronics [6], gravitational lensing [7], smetic C liquid crystals [8], advanced field theory [9], oscillations of buildings in earthquakes [10], etc.In the small-angle regime, the approximation sin θ ≈ θ works and the equation of motion can be linearized, becoming solvable in terms of a sinusoidal function of time.Beyond this regime, the nonlinear nature of the pendulum motion becomes apparent: (i) the period increases with the amplitude; and (ii) the function θ(t) departs more and more from a harmonic behavior.In the nonlinear regime, the simplicity of the sinusoidal solution is lost and the equation of motion becomes unsolvable in terms of elementary functions, so either a direct numerical solution or some analytical approximation is demanded.However, numerical solutions for the pendulum equation loose accuracy for long time intervals due to the accumulation of roundoff errors (inherent to numerical routines) [11], whereas polynomial approximations (including truncated Taylor series) have an intrinsic defect, namely the divergence for very long time intervals, as follows from the fact that lim t→∞ P n (t) = ± ∞ for any polynomial with degree n > 0.Here in this work, I show how these limitations can be circumvented by taking a few terms of the Fourier series expansion of the elliptic function sn(u; k) as the basis to build a simple periodic approximation for θ(t) which is accurate for all t > 0, even for amplitudes near π rad.

Equation of motion and exact solution
A simple pendulum consists of a small bob of mass m suspended by a weightless, thin rigid rod of length fixed at the upper end and free to rotate around the suspension point C, as indicated in Fig. 1.Its oscillatory motion can be studied by applying Newton's second law of motion for rotations, which, in the absence of dissipative forces and taking the counterclockwise direction of motion as positive, reads Figure 1: The simple pendulum motion.At t = 0, it is released at rest from a position that forms an angle 0 < θ0 < π rad with the vertical, then passing by an angle θ, for some 0 < t < T /4, with a scalar velocity v = θ.
where g is the local acceleration of gravity and I = m 2 is the pendulum's moment of inertia.This simplifies to the nonlinear pendulum equation [12,Sec. 4.4]: where ω ≡ g/ .
Beyond this regime, one should retake the nonlinear differential equation in Eq. (2).Alternatively, an expression for the inverse function t(θ) can be derived for any 0 < θ 0 < π rad from the conservation of mechanical energy, without a deeper discussion on differential equations [12][13][14].On taking the zero level of the potential energy at the lowest point, one finds [13] On isolating θ, one has where the '+' sign is for the counterclockwise motion and the '−' sign is for the clockwise one.The integration from θ 0 to any θ(t) with 0 < t ≤ T /2 (hence using the '−' sign) yields This improper integral becomes a proper one making cos β = 1 − 2 sin 2 (β/2), followed by the substitution sin (β/2) = k sin φ, where k ≡ sin(θ 0 /2).This yields where sin ϕ = sin (θ/2)/k, is the complete elliptic integral of the first kind, and is the incomplete one [15].The parameter k is the elliptic modulus and ϕ = arcsin [sin(θ/2)/k] is the amplitude of u(ϕ; k).Though these integrals cannot be expressed in a closed-form in terms of elementary functions only [16], a formula for the exact pendulum period T can be derived by taking t = T /4 and ϕ = 0 in Eq. ( 7), which corresponds to θ(T /4) = 0.This yields The computation of the exact period T then demands the numerical evaluation of K(k), which is usually done via numerical integration or some recursive iterations in a computer [6].Alternatively, one can use some of the practical approximate formulae available in literature (see, e.g., Refs.[17][18][19][20]).Anyway, Eq. (10) shows that the pendulum period increases from T 0 to infinity as the amplitude θ 0 goes from 0 to π, as expected since the top vertical position θ 0 = π is a point of (unstable) equilibrium [14].This increase of T with θ 0 makes the harmonic solution in Eq. ( 3) inappropriate for large-angle oscillations.
Though the exact solution of Eq. ( 2), for any 0 < θ 0 < π rad, cannot be expressed as a finite combination of elementary functions, it can be written in terms of elliptic functions, as proved independently by Abel and Jacobi in 1827 [21,22].For a given k, k 2 < 1, the three basic Jacobi elliptic functions are defined as inverses of the elliptic integral u(ϕ; k) [15]: For real values of u, these functions are smooth and periodic [23].From Eq. ( 9), it promptly follows that du dϕ = 1 dn(u; k) , which allows us to apply the chain rule to get the derivatives of the above functions with respect to u.The results are [11,15]: With these derivatives in hands, it is easy to show that [11] The connection with the pendulum equation arises by making the change of variables z = sin ϕ = sin(θ/2)/k directly in Eq. ( 2), which yields A comparison of Eqs. ( 13) and ( 14) promptly leads to the general solution [14,17]: For the initial conditions adopted in Eq. ( 3), which now read z(0) = 1 and ż(0) = 0, one finds C = 1 and ϕ 0 = K(k), so z(t) = sn (ωt + K(k); k).The exact solution θ(t) then reads: where, for simplicity, K denotes K(k), as will be done hereafter.

Periodic approximate solutions
In order to overcome the disadvantages of numerical solutions and polynomial approximations in describing the pendulum motion, let us restrict ourselves to periodic analytical approximations only.As a first approximation of this kind, I modify the small-angle solution in Eq. (3) to where T is the exact pendulum period, as given in Eq. ( 10).This approximation is slightly better than a naive cosine approximation I did propose a decade ago [17, Eq. ( 16)].There in that work, I did show that such harmonic approximation is practically able from the exact solution θ(t) for amplitudes below π/4 rad, deviating significantly from it only for amplitudes above π/2 rad.This 'loss of harmonicity' has been confirmed by numerical Fourier analysis in a very recent work [24].On searching for better periodic approximations, I have noted that the elliptic function sn(u; k) which composes the exact solution established in Eq. ( 16), is a continuous, periodic function for real values of u, with period 4 K [16,23], so it satisfies the Dirichlet's conditions for convergence, being natural to investigate its Fourier series expansion.On making u = 2 K x/π in the odd function sn(u; k), one gets sn(2 K x/π; k), an odd function with period 2π to which the Fourier sine-series with Fourier coefficients converges for all real values of x.This last integral can be solved analytically using a contour integral in the complex plane, as shown in Ref. [25].The main steps of this procedure are shown in the Appendix.The final result is: where q = e − π K /K and K The fast convergence of this Fourier series suggests that its truncation to a few terms N ≥ 1 can approximate sn(u; k) accurately enough for to be a good approximation for the exact solution θ(t).
In the next section, this possibility will be investigated.
For a not-so-large amplitude such as θ 0 = π/3 rad, the exact solution θ(t) (black, solid curve), the harmonic approximation (red, dotted curve), our Eq.( 17), and the Fourier-based approximation (blue, dashed curve), our Eq.( 22), are plotted in Fig. 2. As may be seen, both approximations work well, the harmonic one (relative error < 1%) being somewhat better than the Fourier-based approximation taken with just one term, whose relative error reaches 2%.However, the harmonic approximation is much poorer than the Fourier one with N = 2 (not included in Fig. 2), whose maximum relative error is only 0.04%.
For θ 0 = π/2 rad, which is the largest amplitude for smooth oscillations of a pendulum sustained by a flexible wire, the results are shown in Fig. 3.As seen in the upper panel, the harmonic approximation θ(t) seems to remain accurate, but the error graph in the lower panel reveals a maximum relative error of 2.3%, much larger than 0.2% obtained with N = 2 in the Fourier-based approximation.
For larger amplitudes, a rigid rod is demanded for producing smooth oscillations and the loss of harmonicity is more striking. 1For instance, for θ 0 = 4 5 π rad, the harmonic approximation is highly inaccurate, whereas the Fourier-based approximation with N = 3 is accurate to a 0.6% level.Even for amplitudes as large as θ 0 = 7 8 π rad (i.e., 157.5 • ), as shown in Fig. 4, the Fourier approximation with N = 4 presents a maximum relative 1 When the mass of the rod is not negligible in comparison to the mass of the bob one has a physical pendulum, for which both the length (i.e., the distance between the center of mass and the suspension point) and the moment of inertia in Eq. ( 1) must be reconsidered.See Ex. 11.2 of Ref. [12].16); the dotted curve (red) is for the harmonic approximation θ(t), Eq. (17); and the dashed curve (blue) is for the Fourier-based approximation θ(t), our Eq.(22), with N = 1 only.Note that the solid and dotted curves are almost indistinguishable, whereas the dashed (blue) curve departs somewhat from that set.Lower panel: deviations from θ(t) with θ(t) (dotted curve) and θ(t) (dashed curve).error of only 0.4%.Even for very large amplitudes such as θ 0 = 0.95 π rad (i.e., 171 • ), the Fourier approximation with N = 6 attains a 0.3% accuracy.

Conclusion
In this note, periodic analytical approximations for the exact solution of the pendulum equation of motion are proposed.As a first approximate solution, I have modified the harmonic approximation proposed in Ref. [17], which resulted in an accurate approximation for amplitudes up to π/3 rad.However, due to the unavoidable loss of harmonicity with the increase of amplitude, the harmonic approximation soon becomes poor.In fact, for amplitudes above π/2 rad, which are of interest for some experiments [4,8,[26][27][28], the deviations are significant, as shown in Fig. 3.Then, I have found it natural to take the periodicity of the exact solution into account to derive an analytical Fourier series expansion, giving continuity to a numerical treatment I have developed (with co-authors) in a very recent work [24].This task revealed many complexities due to the inverse sine function present in the exact solution, our Eq.( 16), so I decided to expand only the (periodic) Jacobi elliptic function sn(u; k), rather than the whole function θ(t).The truncation of this series to only a few terms has revealed itself as a good approximation technique, yielding relative errors smaller than 0.4% even for amplitudes as large as 7 8 π rad with four terms only, as shown in Fig. 4. A previous treatment of this sort is not found in the literature and the advantages of the periodic approximate formulae derived here are certainly their simplicity and accuracy, which makes them useful for analysing large-amplitude oscillations in experiments, working well for arbitrarily long time intervals.This kind of periodic approximations may be of general interest for science and engineering applications in which the natural modes of oscillatory structures must be determined accurately, e.g.avoiding dangerous resonances of buildings in earthquakes [10].In most real oscillatory systems the inclusion of dissipative effects (e.g., friction and air resistance) is necessary for a more realistic modelling, so the extension of the approach put forward here to dissipative pendular systems is left as a suggestion to the more interested readers.

Figure 2 :
Figure 2: Periodic approximations for the exact solution θ(t) of the pendulum equation, for θ0 = π/3 rad.Upper panel: the solid curve (black) is for the exact solution, Eq. (16); the dotted curve (red) is for the harmonic approximation θ(t), Eq. (17); and the dashed curve (blue) is for the Fourier-based approximation θ(t), our Eq.(22), with N = 1 only.Note that the solid and dotted curves are almost indistinguishable, whereas the dashed (blue) curve departs somewhat from that set.Lower panel: deviations from θ(t) with θ(t) (dotted curve) and θ(t) (dashed curve).

Figure 3 :
Figure 3: The same of the previous figure, for θ0 = π/2 rad.The Fourier approximation was taken with N = 2. On the left panel, the solid and dashed curves are practically indistinguishable, but the dotted (red) curve already departs from them.

Lima e20180202- 5 Figure 4 :
Figure 4: The same of the previous figures, for θ0 = 7 8 π rad.The Fourier approximation was taken with N = 4. Again, the solid and dashed curves are practically indistinguishable, whereas the dotted (red) curve departs much from them.