Fourier series expansion for nonlinear Hamiltonian oscillators

The problem of nonlinear Hamiltonian oscillators is one of the classical questions in physics. When an analytic solution is not possible, one can resort to obtaining a numerical solution or using perturbation theory around the linear problem. We apply the Fourier series expansion to find approximate solutions to the oscillator position as a function of time as well as the period-amplitude relationship. We compare our results with other recent approaches such as variational methods or heuristic approximations, in particular the Ren-He’s method. Based on its application to the Duffing oscillator, the nonlinear pendulum and the eardrum equation, it is shown that the Fourier series expansion method is the most accurate.


I. INTRODUCTION
The dynamic behavior of nonlinear oscillators has been traditionally analyzed by using perturbation methods ͑such as the Lindstedt-Poincaré method or the Shohat expansion͒ where the initial amplitude of oscillation is assumed to be a small parameter, or considering the averaging method ͓1͔.More recently there have been an amount of works devoted to improve the previous methods with alternative analytical techniques that do not need to define a small parameter, such as variational principles ͓2-4͔ or variational iteration methods ͓5͔.The variational principle is constructed to obtain the period-amplitude relationship by maximizing a functional over a set of trial functions.The method ensures the existence of a trial function, which yields the exact result but does not provide it.The most common approach is to make use of the trial function that gives the exact result for the linear problem.The variational iteration method consists in constructing a correction functional by a Lagrange multiplier, which can be identified via the variational theory.The method also uses the solution to the linear problem as initial approximation and, although we have not been able to prove it generally, this is the reason why both variational methods seem to give the same results for the period-amplitude relationship.Other more recent methods are the homotopy perturbation ͓6͔, the harmonic balance method ͓7͔, and heuristic methods based on ancient Chinese mathematics ͑namely, Ren-He's method͒ ͓8͔.
Fourier series expansion is a useful tool for solving linear differential equations.However, this method is also adequate to find approximate solutions to nonlinear problems such as population pattern formation in continuum media ͓9,10͔ or to characterize chaotic populations in lattices living in finite domains under Dirichlet boundary conditions ͓11͔.
In this work, we apply the Fourier series expansion method to solve nonlinear Hamiltonian oscillators and compare our results with numerical calculations and with the approximated solutions supplied by some of the aforementioned methods.We rigorously prove that the leading term of the expansion recovers the solutions obtained from the variational principle developed in ͓2͔ and, as it occurs in the virial expansion for a real gas, each of the terms in the expansion is of lower order than the previous one.Hence, we stress that Fourier series expansion goes beyond the linear approach on which the variational principles are based.In consequence our method is in better agreement with the numerical solutions.In some sense, the Ren-He's method also makes use of the solution to the linear oscillator to get a new solution in a kind of iterative procedure.In fact, Ren-He's method is nothing but the first iterate of the classical Picard iteration method.We also observe that the Ren-He's method shows pathological results in some family of nonlinear oscillators due to the appearance of secular terms that diverge with time.In Secs.II and III, we discuss these general results and in Sec.IV, we provide comparison with three widely known examples: the Duffing oscillator, the nonlinear pendulum and the eardrum equations.

Let us consider nonlinear Hamiltonian oscillators of the form
where a is the initial amplitude of motion.The prime symbol stands for the temporal derivative, is a parameter and f͑u͒ is a nonlinear function of the position u.When = 0 the oscillator is linear, u͑t͒ = a cos͑t͒ and T =2.The Fourier series method consists in proposing a solution to Eq. ͑1͒ of the form u͑t͒ = ͚ m ␣ m m ͑t͒, where the coefficients ␣ m are to be determined and m ͑t͒ is a basis of orthogonal functions which must satisfy the initial conditions Eq. ͑2͒.To this end, we choose the trigonometric basis m ͑t͒ = cos͑k m t͒.It is clear that uЈ͑0͒ = 0 for any k m ; however, from the definition of period T, one has u͑T / 4͒ = 0, so that it is necessary that k m T / 4=͑2m −1͒ / 2 for m =1,2,3,....In consequence, we propose To find the solution we need to obtain the set of coefficients ␣ m .Because of the orthogonality of the trigonometric func-tions, one gets a set of nonlinear algebraic equations for the coefficients ␣ m by introducing Eq. ͑3͒ into Eq.͑1͒ and then multiplying this equation by cos͓2͑2n −1͒t / T͔ and integrating over t from 0 to T / 4. As this system of equations cannot be solved analytically, we consider the following approximation: let us assume that in the expansion Eq. ͑3͒ each term is smaller than the previous one in such a way that the main contribution to the solution Eq. ͑1͒ comes from the first term, that is ␣ 1 cos͑2t / T͒.Then, we consider f͑u͒ Ӎ f͓␣ 1 cos͑2t / T͔͒, neglecting higher order terms.The set of algebraic equation reduces to cos͓͑2n − 1͒s͔f͓␣ 1 cos͑s͔͒ds.

͑4͒
On setting n = 1 in Eq. ͑4͒ we get the equation to solve ␣ 1 : which can be introduced into Eq.͑4͒ to find ␣ n as a function of the period T.However, what we aim to find is how u͑t͒ depends on the initial amplitude a, so that we need to obtain first the period-amplitude relationship.From Eqs. ͑2͒ and ͑3͒ this relation reads where ␣ 1 is given by Eq. ͑5͒.To prove the validity of the approximation f͑u͒Ӎ f͓␣ 1 cos͑2t / T͔͒ the first term in the expansion Eq. ͑3͒ must be the leading term, and in consequence the condition lim ␣ 1 →0 ␣ n / ␣ 1 = 0 must be satisfied.This is easy to check from Eq. ͑4͒.As f͑ • ͒ is a nonlinear function, when ␣ 1 → 0 one has ␣ n ϳ f͑ • ͒ϳo͑␣ 1 ͒ and lim ␣ 1 →0 o͑␣ 1 ͒ / ␣ 1 = 0.The approximation is more accurate when the amplitude a of motion is small and, in consequence, the linear term becomes more important than nonlinear terms.
Let us now compute the main contribution in the truncated expansion Eq. ͑3͒ up to the first order.In this case The integral in Eq. ͑8͒ a can be written as dz under the change of variable z = a cos͑s͒.Finally, the periodamplitude relationship is given by This result, together with Eq. ͑7͒, proves that the Fourier series expansion up to the first order coincides with the predictions from the variational principle ͓2͔.

III. REN-HE'S METHOD
The Ren-He's method ͓8͔ consists in choosing as the trial function u 0 = a cos͑2t / T͒, and introducing it into Eq.͑1͒ to get uЉ + u 0 + f͑u 0 ͒ = 0, which is a linear ordinary differential equation that can be solved by integrating twice under the initial conditions Eq. ͑2͒.The solution reads The period-amplitude relationship is obtained from Eq. ͑10͒, through the condition u͑T / 4͒ = 0, to yield the implicit relationship This is a heuristic method with many known variations that fits well with the exact relation between period and amplitude.However, less attention is paid in checking this method by comparing the result Eq. ͑10͒ with the numerical or exact solution.We now show that this method fails when applied to nonlinear oscillators of the form uЉ + u + u n = 0 with n an even number.Inserting f͑u͒ = u n into Eq.͑10͒ one gets The double integral in Eq. ͑12͒ can be explicitly calculated.
When n is an even number while for odd n one finds ͓13͔ ͵ As can be seen, the main difference is that when n is even a term proportional to t 2 appears, a secular term that cannot be removed.This pathological behavior leads u͑t͒ to decrease toward −ϱ as t → ϱ.As we will see in an example, the agreement between the solution of the Ren-He's method and the numerical solution is bad, even for relatively short times.

IV. EXAMPLES
In this section, we apply the Fourier series, variational principle and Ren-He's methods to find the solution of Eq. ͑1͒ and the period-amplitude relationship for different functions f͑ • ͒.The results obtained for the oscillator position from the three approaches will be compared to numerical solutions performed by using the Runge-Kutta method RK45.The period-amplitude relationship can be also obtained by integrating Eq. ͑1͒.The velocity of the oscillator is defined by v =−uЈ in such a way that v Ͼ 0 for t ͓0,T / 4͔.After a first integration and taking into account the initial condition v͑0͒ = 0 one finds . ͑15͒ Integrating Eq. ͑15͒ and using u͑T / 4͒ = 0 the periodamplitude relationship is For many nonlinear oscillators the period given by Eq. ͑16͒ can be expressed in terms of elliptic integrals, hypergeometric functions, or other special functions ͓1͔.However, we will solve Eq. ͑16͒ numerically in each example to check the accuracy of the methods here discussed.

A. Duffing oscillator
The well-known Duffing oscillator is described by Eq. ͑1͒ with f͑u͒ = u 3 .From Eq. ͑4͒ we have where ␦ nm is the Kronecker delta.The period-amplitude relationship Eq. ͑6͒ is expressed as The solution for the position of the Duffing oscillator is, from Eqs. ͑3͒ and ͑17͒, u͑t͒ = 2 where T is given by Eq. ͑18͒.
The solution obtained from the variational principle is simpler: from Eq. ͑9͒, the period-amplitude relationship is given by which coincides with the result found in ͓2͔.From Eqs. ͑7͒ and ͑20͒, the solution to the position is

͑21͒
Now let us write the solutions obtained from the Ren-He's method.From Eqs. ͑10͒ and ͑11͒ one gets respectively.Both results have been recently obtained by Ren and He ͓8͔.
In Fig. 1, we plot the position of the Duffing oscillator as a function of time.As shown in Fig. 1͑a͒, the three methods agree very well in absolute terms.Figure 1͑b͒ is an amplification of the same Fig.1͑a͒, and allows us to see that the Fourier series expansion result ͑squares͒ is in better agreement with numerical results ͑solid lines͒ than variational ͑dashed lines͒ or Ren-He's ͑dotted lines͒ methods.The difference between the approximate methods and the numerical solution grows with time, as one can appreciate in Fig. 1͑b͒.
In Fig. 2, we plot the period-amplitude relationship for the Duffing oscillator.Similarly, in Fig. 2͑a͒, we appreciate that all three methods agree very well with the numerical results obtained ͑solid line͒ in general terms but if we look in deeper detail as done in Fig. 2͑b͒, it is found that the Fourier series method is again the best fit to the numerical result, followed by the variational approach.This means that Fourier series expansion slightly improves the first order ͑leading͒ contribution obtained from the variational principle.

B. Nonlinear pendulum
The nonlinear Eq. ͑1͒ for the pendulum reads uЉ + sin͑u͒ = 0, where u͑t͒ describes the angle as a function of time so that the maximum amplitude of motion is ͑hence, a Ͻ ͒.Although this case is somewhat different from model Eq.͑1͒, it is a very instructive example.The Fourier series expansion yields where J n ͑ • ͒ is the Bessel function of first kind of order n.
Taking n = 1 we obtain the transcendental equation 2 2 ␣ 1 = T 2 J 1 ͑␣ 1 ͒ for ␣ 1 , which allows us to find ␣ 1 = ␣ 1 ͑T͒.Introducing this result into Eq.͑24͒ we get ␣ n = ␣ n ͑T͒ for n Ͼ 1.The period-amplitude equation is then acquired by summing all ␣ n , hence FIG. 1. Position of the Duffing oscillation as a function of time for = a = 1.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.In panel ͑a͒ time ranges from 0 to 20 and in panel ͑b͒ it ranges over a small period of time in order to appreciate the differences between the different approaches.The time is given in dimensionless units.FIG. 2. plot for the Duffing oscillator =1.T and a are given in dimensionless units.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.In panel ͑a͒ we show the range 0 Ͻ a Ͻ 50 while in panel ͑b͒ we show an enlargement to appreciate the differences between methods.
and the solution for the angle as a function of time reads

͑26͒
Truncating the expansion at first order, and since as ␣ 1 Ӎ a, one has the result for the period-amplitude relationship obtained from the variational principle ͓4͔ and the solution for the angle

͑28͒
The equation to solve for the Ren-He's method is uЉ + sin͑u 0 ͒ = 0, under the initial condition Eq. ͑2͒ has the solution

͑29͒
where we have considered the property Setting u͑T / 4͒ = 0 into Eq.͑29͒ we find the formula

͑30͒
for the period-amplitude relationship.
As in the Duffing oscillator, all methods agree quite well with the numerical solutions for the angle as function of time ͓see Fig. 3͑a͔͒, once again the Fourier series expansion provides the best fit, followed by the variational expansion ͓see Fig. 3͑b͔͒.When the initial amplitude a tends to its maximum value the period should tend to infinity.This can be observed in Fig. 4͑a͒.As a → − all methods tend to infinity but their fit to the numerical solutions is worst accurate.For a Ͻ / 2 all methods are in good agreement and, as we can see in Fig. 4͑b͒, the Fourier series expansion gives the best results, followed by the variational result.However, it diverges faster than the other two methods.This is due to the presence of a separatrix at u = ͑see Appendix͒.

C. Eardrum equation
The equation of motion for the human eardrum can be written by Eq. ͑1͒ with f͑u͒ = u 2 ͓1͔.The Fourier series expansion predicts, from Eq. ͑6͒, the period-amplitude relationship The position of the oscillator is, from Eq. ͑3͒, given by FIG. 3. Angle of the pendulum as function of time for a =1.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.The time is given in dimensionless units.

͑32͒
Truncating the Fourier series expansion up to first order and considering ␣ 1 Ӎ a we recover the result of the variational principle.The period-amplitude relationship reads from Eq. ͑9͒ and the position is, from Eq. ͑7͒, u͑t͒ = a cosͩtͱ1 + 8a 3 ͪ.

͑34͒
The period-amplitude relationship obtained with the Ren-He's method ͓Eqs.͑10͒ and ͑11͔͒ is and the position Due to the presence of a separatrix and heteroclinic orbits for a Ͼ 1 / 2 ͑see Appendix͒ the methods are able to fit the numerical results only for a Ͻ 1 / 2. In Fig. 5͑a͒ it is shown that the Ren-He's method fails to agree with the numerical results for t Ͼ T due to the secular term in Eq. ͑36͒, while the other two methods agree quite well.In the case depicted in Fig. 5 the contribution of the terms with n Ͼ 1 in the Fourier series expansion Eq. ͑32͒ does not improve substantially the firstorder correction as we can see from Fig. 5͑b͒, probably due to the fact that the initial amplitude is small ͑in this figure a = 0.1͒.The Fourier series expansion and the variational method give almost the same results.In Fig. 6, it is shown that the Ren-He's method provides a reasonably good method for the period-amplitude relationship in contrast to its prediction for the oscillator position.Nevertheless, this time the best method is by variational principle.

V. CONCLUSIONS
We have analyzed different methods for solving the problem of nonlinear Hamiltonian oscillators, emphasizing the Fourier series expansion as an accurate method, where an appropriate basis of orthogonal functions, that accomplish the boundary conditions, permits to write the oscillator position as a sum of different terms.The comparison of the position as a function of time and of the period-amplitude relationship is better than with variational approaches and the recent Ren-He's method.Therefore, the Fourier series expansion generally improves the other results, as it has been shown in the different examples studied: Duffing oscillator, nonlinear pendulum, and eardrum equation.In fact, we have FIG.4. Period-amplitude plot for the nonlinear pendulum for = 1.In panel ͑a͒ we show the complete range of the initial amplitude ͑0 Ͻ a Ͻ ͒.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.In panel ͑b͒ we show only small values for the initial amplitude.T and a are given in dimensionless units.
shown that the Fourier series expansion is a higher order correction to the variational approach, which is recovered truncating the Fourier expansion at first order, and all the other terms approximate better the real result, as in a kind of virial expansion.A strong limitation for all analytical methods is the presence, for some oscillators, of a separatrix, ͑i.e., of heteroclinic orbits͒.The Ren-He's method is simple and convenient only in some situations, as it fails in predicting the position of nonlinear oscillators with nonlinear terms of powers of even exponents although the system is far from the separatrix.The Fourier series expansion is rather general, and specially accurate for small amplitudes.Although this method is shown to be more accurate than the other methods, for the case of the nonlinear pendulum it diverges faster near the separatrix.It would be interesting to study if this method holds for dissipative oscillators but it needs a more careful analysis and could be the subject of future works.

APPENDIX: SEPARATRICES
Equation ͑1͒ can be integrated and one gets

͑A1͒
where E is an integration constant.Since Eq. ͑1͒ is conservative, E can be regarded as the total energy, i.e., the sum of the kinetic energy K = uЈ 2 / 2 and the potential energy

͑A2͒
The fixed points of the system, u 0 , satisfy the equation u 0 + f͑u 0 ͒ = 0, i.e., they are solutions to ͑dV / du͒ u=u 0 = 0. Their FIG. 5. Position of the eardrum oscillator as function of time for = a = 0.1.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.The time is given in dimensionless units.FIG. 6. Period-amplitude plot for the eardrum with =1.T and a are given in dimensionless units.The numerical solution ͑solid curves͒ is depicted together with the other three methods studied: Fourier series expansion ͑squares͒, variational principle ͑dashed curves͒, and the Ren-He's method ͑dots͒.In panel ͑a͒ we show the range 0 Ͻa Ͻ 1 / 2, for which the oscillator describes homoclinic orbits.In panel ͑b͒ we show an enlargement to appreciate the differences between methods.stability is given by the sign of the second derivative of the potential.If ͑d 2 V / du 2 ͒ u=u 0 Ͼ 0 then u 0 is stable and if ͑d 2 V / du 2 ͒ u=u 0 Ͻ 0 then u 0 is unstable.If the system Eq.͑1͒ has only one fixed point for which V is minimal ͑as occurs for the Duffing oscillator͒ it is a center and its dynamics is always described by homoclinic orbits ͑each orbit has a different value of the energy E͒, that is, the system oscillates around the fixed point.However, the nonlinear pendulum and the eardrum equation have a center and a saddle point.For example, consider the nonlinear pendulum equation with −2 Յ u Յ 2. u 0 = 0 is a center while u 0 = and u 0 =− are saddle points.This means that there exists a region called separatrix that is a frontier between the attraction basins of both fixed points or, in other words, is a frontier between homoclinic and heteroclinic orbits.For the case of the ear-drum equation there is a center at u 0 = 0 and a saddle point at u 0 =−1/ .Now the question is to find where the separatrix is located.In particular, we are interested in finding the values for the initial amplitude a for which the system exhibits both homoclinic or heteroclinic orbits.Let E ‫ء‬ be the energy level corresponding to the saddle point and V min the energy level corresponding to the center.If V min Ͻ E Ͻ E ‫ء‬ , the orbits are homoclinic and if E Ͼ E ‫ء‬ they are heteroclinic.Hence, the separatrix corresponds to E = E ‫ء‬ .For the pendulum f͑u͒ = ͓sin͑u͒ − u͔ / and from Eqs. ͑2͒ and ͑A1͒ one finds E = −cos͑a͒ and E ‫ء‬ = 1, so that the orbits are homoclinic for 0 Ͻ a Ͻ .For the eardrum equation one has E = a 2 / 2+a 3 / 3 and E ‫ء‬ =1/ 6 2 .In this case, the orbits are homoclinic for a Ͻ 1 / 2.