Improvement of Harmonic Balance Using Jacobian Elliptic Functions

We propose a method for finding approximate analytic solutions to autonomous single degree-offreedom nonlinear oscillator equations. It consists of the harmonic balance with linearization in which Jacobian elliptic functions are used instead of circular trigonometric functions. We show that a simple change of independent variable followed by a careful choice of the form of anharmonic solution enable to obtain highly accurate approximate solutions. In particular our examples show that the proposed method is as easy to use as existing harmonic balance based methods and yet provides substantially greater accuracy.


Introduction
Ordinary differential equations (ODEs in short) are ubiquitous in fundamental science as well as in engineering. Indeed they commonly arise as direct results of the application of some fundamental laws in the various fields of science or engineering. A classical example here is Newton's second law of motion. They also often arise indirectly, for example, in the intermediate steps of solving other types of problems such as partial differential equations (PDEs). Solving ODEs, especially analytically, thus appears to be of great importance for gaining insights into real-world or engineering problems. This is yet a very challenging task since the ODEs of interest, being usually nonlinear, are rarely susceptible to exact analytical solutions. In fact, the lack of a general and systematic methodology for solving these nonlinear equations is probably the most important difficulty in the determination of their analytical solutions. So, for the important class of ODEs associated with oscillatory behaviors, numerous techniques which enable to obtain some approximations to the desired solutions have been developed. These techniques may be classified in two groups: perturbation and non-perturbation.
The most famous of perturbation approaches include the Lindstedt-Poincaré (LP) method, the method of multiple scales (MMS) and the Krylov-Bogoliubov-Mitropolsky (KBM) method. These now classical methods are known to have some serious limitations. For instance they are useless for equations which describe essentially nonlinear oscillators, such as (1) below with In contrast, non-perturbation techniques of which the method of harmonic balance (HB) is a well-known example, do not suffer these limitations. But the straightforward application of this last method leads to systems of nonlinear algebraic equations for the coefficients of the truncated Fourier series assumed for the solutions; which are still very difficult to solve. However the acuteness of this difficulty can now be reduced considerably thanks to the idea proposed initially by Wu and Li [1] and further refined by Wu and co-workers [2] [3]. It consists of linearizing the governing equations before applying the HB itself, which then leads to linear algebraic equations instead of nonlinear ones. The resulting method, harmonic balance with linearization (HBwL), has been applied to various types of conservative symmetric oscillators and appears to be both quite efficient and very simple [2] [3], and references therein. Considering these two interesting properties, the method of HBwL is a good candidate for improvements which can make it to be more versatile than just conservative or symmetric systems. In this respect it has been adapted to handling dissipation terms in ODEs [4], and more recently to solving asymmetric oscillator equations [5].
Another important class of non-perturbative techniques involves the approximation of the nonlinear restoring force in a given oscillator ODE by some simple forms for which the exact solution can be readily obtained. The most common technique in this class is the method of equivalent linearization. In this approach the approximate ODE has the simple form of a linear harmonic oscillator equation [6]. Besides the method of equivalent linearization is the so-called method of "cubication" in which the restoring force function is expanded up to third order in Chebyshev polynomials. The method of cubiccation thus approximates a given ODE by a cubic Duffing equation whose exact solution is the approximate solution to the original ODE. It is well known that the exact solution of the Duffing equation can be expressed in explicit form using Jacobian elliptic functions [7] [8].
In this paper, we investigate a further generalization of the HBwL by using Jacobian elliptic functions instead of circular trigonometric functions [7]. Our objective is to develop a higher order approximation than the single cnoidale elliptic function which is obtained when a given ODE is first approximated by a Duffing equation. We achieve this goal through a simple change of variable based on the properties of the Jacobian elliptic functions, which results to an ODE in which the restoring force is not approximated and which is solved approximately following the method of HBwL or other generalizations of the method of HB such as the rational HB.
The paper is organized as follows. The main idea sustaining the elliptic harmonic balance is introduced in the next section. The illustrations are presented next in Section 3, where examples of oscillators with polynomial and rational restoring force functions are considered. Section 4 contains the conclusion of the paper.

Presentation of the Approach
We consider a nonlinear second-order differential equation of the general form which governs the time (t) evolution of the state of the state (x) of a conservative single degree-of-freedom system. Here, an overdot denotes differentiation with respect to time t. Thus,

( )
x x t = is a real function. We assume without loss of generality that 0 A > . The restoring force function f is a general nonlinear function of x which is required to be continuous and to satisfy ( ) ; so that the system oscillates symmetrically around 0 x = . Our objective is to determine an approximate analytic solution to the periodic solution of the above equation. It is customary in classical techniques [3] [9]- [12] to introduce a time scaling according to 2π , .
Here, T and w are the unknown time-period and angular frequency of the sought periodic solution. Under this scaling, (1) becomes where we have set 2 w Ω = , and a prime denotes a differentiation with respect to τ. The periodic solutions of this latter equation are 2π-periodic in τ.
In this paper we use a different time-scaling than (2) which is simply linear in τ. Following Yuste Bravo [7], we consider the generally nonlinear transformation given by where am is the Jacobian elliptic amplitude function of parameter m while K is the complete elliptic integral of the first kind. In this paper we adopt the notation of [13] for designating elliptic functions. The parameter m is an additional unknown to be determined as part of the solution. It should be noted that (2) is a particular case of (4), corresponding to the fixed value m = 0. Using the following well-known property of the am function one can easily show that the nonlinear time-scaling (4) transforms (1) Notice that (6) reduces to (3) for m = 0. The former differs from the latter only by two terms which are parametrically linear in the derivatives. Thus if the HB can be applied to the latter, which requires that one be able to compute the Fourier series of ( ) for a given truncation of the 2π-periodic Fourier expansion of ( ) x τ , then it should equally be applicable to (6). Once the solution is obtained in terms of τ, its expression is transformed into a polynomial in  x Then, by making the substitutions sn wt m τ → one can express the approximate solution in terms of t. Here cn and sn are respectively the Jacobian elliptic cosine and the Jacobian elliptic sine functions [13] [14]. It is worth to mention that in basic introduction to Jacobian elliptic functions it is supposed that 0 1 m ≤ ≤ . However this assumption is not a real restriction for our investigation. In effect, as the approximate solutions considered here are ultimately expressed in terms of cn and sn, the following transformations ( [13], p. 573) ( ) are applicable respectively for m < 0 and m > 1. Thus the results that would be obtained using the proposed approach are valid for all real values of such that 1 m ≠ (for this valueof , the Jacobian elliptic functions become hyperbolic and nonperiodic). Moreover, the Jacobi's imaginary transformation ( [13], p. 574) in which 2 1 i = − can be applied when the value of 2 w Ω = is negative; as long as the expression of the approximate analytical solution to the ODE involves only cosine terms. Note that in (8)-(10),

Illustrations
Two examples are now presented for illustration.

The Cubic-Quintic Duffing Oscillator
Our first example involve the cubic-quintic Duffing oscillator described by (1) with f defined such that , the oscillations will be symmetric within the bounds [ ] , A A − . It is well known that the Fourier expansion of ( ) x τ contains only odd harmonics in this case [15]. The simplest truncation of this expansion which satisfies the prescribed initial conditions consists of the fundamental harmonic only, Substituting this in (6) 16 15 The explicit solution associated with (12) is as in (18a)  The results obtained at this level are exactly the same that would be obtained by cubication [16]. In effect, the third order Chebyshev expansion of (11)  16 where 1 T and 3 T are the Chebyshev polynomials of degrees 1 and 3 respectively: The method of cubication therefore approximates (1) with (11) by The solution of (17a) can be written explicitly as [16] ( ) ( ) ; x t Acn wt m One can easily verify that for μ and λ defined by (17b) the expressions of and in (18b) and (14) are the same. This shows that the method of cubication and the lowest order harmonic balance in our proposed approach are equivalent. But unlike the method of cubication whose results would be limited to the above, the transformation proposed herein allows for further improvments. For instance, the approximate solution to (6) can be sought directly in the form The required Fourier series expansion is quite easily calculated for polynomial restoring force functions. By equating the coefficients of ( ) The coefficients of the above polynomial equations are as follows: 4 as was suggested by Yuste Bravo [7]. The expression of b following the linearization procedure is found in this case to be ( ) As can be seen, this expression is independent of the coefficient of the linear term α in contrast to (23). The consequence of this is that (25) is less versatile than (19) to accurately approximate the solutions of the Duffing oscillator for all combinations of the coefficients of its terms. In effect, we have noticed that (25) is very poor when 0 β = as can be observed in Figure 1. Even worse, it becomes completely divergent for some combinations of parameters values and oscillation amplitude A. This behavior can be understood by observing that the Fourier series of the residual of (6) with ( ) ( ) and Ω and m given by (14) does not contain the third harmonic. We may conclude that in the HBwL in general, the harmonics to include in the correction to a given stage should be higher than or equal to the least harmonic of the residual terms of that stage. Thus in the remaining part of this work, we consider only the most versatile expression (19) in our discussion.
To further appreciate the accuracy of the various approximations obtained above, we use the relative error Calculations using the approximate expressions derived above show that the relative errors of the first-order and second-order analytical approximations as compared to the exact solution are respectively less than 0.37% and 0.072%. In comparison it would be necessary to proceed up to the third-order approximation of the standard HBwL to obtain a result just better than our least accurate first-order result; with a relative error of 0.23% [17].
In Figure 2 we compare the relative errors of the single harmonic approximation or cubication, the standard HBwL, and the anharmonic solution when 1 α β δ = = = . All the three approximations are highly accurate for small values but yet of order O(1) of A. For larger values of the amplitude, the second order approximation of our present approach appears to yield the smallest relative error of 0.072% for the whole range of allowed oscillation amplitudes for the parameter values chosen. This last observation concerning large is the same for all combination of ( ) , α β with 0 δ > .
A further comparison is presented in Figure 3 for situations where heteroclinic orbits connecting two distinct equilibrium points exist. Once again we observe that our approach generally provides the best accurate approximation. We should point out however that as the oscillation amplitude approaches critically the value for which heteroclinic trajectory is realized, the accuracy of the standard HBwL supersedes our result.

An Oscillator with Rational Restoring Force
As a second example, we consider a conservative nonlinear oscillatory system in which the restoring force has a rational form. Specifically we choose ( ) 2 . 1 A peculiarity of this function which is worth noting is that it is not dominated by odd-power monomials in both of the limiting cases 0 A → and A → ∞ in contrast to the cubic-quintic Duffing oscillator or other rational restroring force oscillators such as the Duffing-harmonic oscillator [16]. In effect while ( ) f x x ≈ for small amplitude oscillations, on has for large amplitude oscillations.
It follows that the expressions of the parameters m and w of the approximate solution of the problem using the method of cubication are 4 To appreciate the accuracy of the two approximate results in (33) and (34)  as a function of A. It appears clearly from this figure that the result obtained from the method of cubication so diverges from the exact result for large oscillation amplitudes that it can simply be considered invalid. However a good match is observed between the result of our proposed approach and the exact result. Figure 5 indicates that the maximum relative error on the period is achieved for A → ∞ and is less than 1.6%.

Conclusion and Remarks
In this paper we have investigated the approximation of periodic solutions to autonomous single degree-offreedom oscillators equations using the Jacobian elliptic function with the objective of improving the method of cubication. To this end we have shown that the properties of these functions can be exploited to put such ODEs  in a form for which the standard HB and like methods can be applied, while the result is however expressed in terms of Jacobian elliptic functions. Our investigation reveals that in the HBwL the harmonics to include in the correction to a given stage should be higher than or equal to the least harmonic of the residual terms of that stage.
With the change of variable sustaining our approach, the analysis can be carried out to higher order, just as in the standard HBwL. For our examples, comparison to the standard harmonic balance indicates that our approach is generally the best, except for periodic orbits too close to a heteroclinic orbit. However, even in this case, the standard HBwL has to be carried to higher order than our approach to gain this advantage. Such a higher order analysis can be quite intricate for non polynomial restoring force. A further interesting property of the solution resulting from our approach is that it also approximates all the harmonics of the exact solution. While the method of rational harmonic balance and the method of cubication also offer solutions with this property, they do not include a mean to carrying the approximations to higher orders as in the present work. In fact our method encompasses the method of cubication which is equivalent to ist first order application.