_{1}

Duffing equation,

There are numerous scientific and engineering phenomena encountering oscillations. The simplest oscillations are described by a generic linear differential equation, ξ ¨ ( t ) + ξ ( t ) = 0 , where ξ(t) is the time-dependent quantity of interest. For instance, ξ(t) could describe the oscillation of the electric charge, q(t), in a LC series electric circuit, or it may describe oscillating coordinate, x(t), of a particle under the influence of a linear force. This equation modifies if one includes retarding, ξ ˙ ( t ) , term and external forces, f(t), resulting, ξ ¨ ( t ) + ξ ˙ ( t ) + ξ ( t ) = f ( t ) , preserving its linearity. Motivation of considering Duffing type equation stems from the fact that the nonlinearity implicitly is embedded in the former equation. For instance, for a mechanical system one might encounter an equation such as, θ ¨ ( t ) + sin [ θ ( t ) ] = 0 . Replacing the second term with its approximate form

yields, θ ¨ ( t ) + θ ( t ) − 1 6 θ 3 ( t ) ≃ 0 . Meaning Duffing equation is implicit to

oscillating motion. As in another instance a nonlinear force may result an equation of motion, x ¨ ( t ) + x ( t ) + ϵ x 3 ( t ) = 0 . Solving these ODEs is a prime interest. Aside searching for their exact solutions [

It is one of the objectives of our investigation to examine the validity of the method [

With these motivations we craft this note; it is composed of four sections. Aside the Introduction, in Section 2 first we brief on the exact and Power Series solution of Duffing equation. Section 3 is the results. Section 4 is the conclusions embodying suggestions for extensions and modifications.

To establish the basis, we begin with the standard Duffing equation,

x ¨ ( t ) + x ( t ) + ϵ x 3 ( t ) = 0 , (1)

Its exact solution with two initial conditions, x ( 0 ) = a 1 and x ˙ ( 0 ) = 0 , is [

x ( t ) = a 1 c n [ u , k ] , (2)

where, u = a 2 t , with a_{2} and k subject to,

a 2 = 1 + a 1 2 ϵ and k = a 1 2 ϵ 2 ( 1 + a 1 2 ϵ ) , (3)

In Equation (2), c n [ ζ , m ] , is the Jacobi elliptic function, c n [ ζ , m ] = cos φ , with φ = F − 1 [ F [ φ , m ] , m ] , where, F [ φ , m ] is subject to [

F [ φ , m ] = ∫ 0 φ d ϑ 1 − m 2 sin 2 ϑ , (4)

As one expects Duffing equation describes oscillations. The severity of deviation of the oscillations given by Duffing equation from the corresponding linear case (ϵ = 0) depends on the value of ϵ. The value of ϵ implicitly impacts the period

of the oscillations given by, T = 4 a 2 F [ π 2 , k ] .

Having established the fundamentals now we cross-examine the accuracy of the Power Series approximation Method outlined in [

The fundamental underlying concept of the approximation method proposed in [

As is outlined in [

x ( t ) = ∑ n = 1 N a [ n ] [ sin ( ω t ) ] n − 1 , (5)

A change of variable, τ = sin ( ω t ) , replaces t, t ≥ 0 , with − 1 ≤ τ ≤ + 1 yielding Equation (1),

ω 2 ( 1 − τ 2 ) d 2 d τ 2 x ( τ ) − ω 2 τ d d τ x ( τ ) + x ( τ ) + ϵ x 3 ( τ ) = 0 , (6)

Noticing significant changes; e.g. Equation (6) gained a dissipative term with a τ-dependent coefficient.

Realizing, x ( τ ) = ∑ n = 1 N a [ n ] τ n − 1 and x 3 ( τ ) = ∑ n = 1 N b [ n ] τ n − 1 , Equation (6) yields the recurrence relationship between the coefficients,

a [ n + 2 ] = a [ n ] [ ω 2 ( n − 1 ) 2 − 1 ] − ϵ b [ n ] ω 2 n ( n + 1 ) , for n = 1 , 2 , 3 , ⋯ (7)

For the chosen initial conditions, x [ t = 0 ] = A = a [ 1 ] and x ˙ [ t = 0 ] = 0 = a [ 2 ] , Equation (7) gives the remaining coefficients yielding the solution of Equation (5). Because of the second initial condition all the even expansion coefficients are zero limiting the polynomial to even powers of the base only. More on this in Conclusion section.

The lack of dissipative term in Duffing equation signifies the conservation of energy. The energy of the system, E = T + V , stays constant. The T and V for a point-like particle of mass m under the action of a nonlinear force, F = k ϵ x 3 , respectively are,

T = 1 2 m x ˙ 2 ( t ) and V = 1 2 k x 2 ( t ) + 1 4 k ϵ x 4 ( t ) , (8)

Equation (8) for the generic case, m = k = 1 is associated with the Lagrangian, L ( x , x ˙ ) = T − V . Applying Euler-Lagrange equation, d d t ∂ L ∂ x ˙ − ∂ L ∂ x = 0 , yields the equation of motion, Equation (1).

Utilizing the time independence of energy and implied initial conditions yields an equation conducive the needed, ω, in Equation (5).

The maximum potential energy occurs at t=0 when the particle is at maximum distance from the equilibrium, i.e. V max ( t = 0 ) = 1 2 A 2 + 1 4 ϵ A 4 . The maximum kinetic energy occurs when the particle passes the equilibrium, T max = 1 2 ω 2 ( 1 − τ 2 ) [ d d τ x ( τ ) ] 2 , at τ = ± 1 2 . Equating these two gives the ω.

T max = V max , (9)

With these outlines we craft a Mathematica program leading the quantity of interest. For objectively chosen numeric values for initial displacement, A, stiffness, ϵ, we investigate the accuracy of the Power Series Method vs. the exact solution; these are given in Sect. 3.

For the chosen value of A and ϵ the number of terms, N, in Equation (5) ought to be determined warrants the convergence of series. To compare our CAS approach to [

30 ω 2 a [ 7 ] = ( − 1 + 16 ω 2 ) ( − 3 A 2 ϵ ( − A − A 3 ϵ ) 2 ω 2 + ( − A − A 3 ϵ ) ( − 1 + 4 ω 2 ) 2 ω 2 ) 12 ω 2 − ϵ ( 3 A ( − A − A 3 ϵ ) 2 4 ω 4 + A 2 ( − 3 A 2 ϵ ( − A − A 3 ϵ ) 2 ω 2 + ( − A − A 3 ϵ ) ( − 1 + 4 ω 2 ) 2 ω 2 ) 4 ω 2 ) (10)

Evaluation of the kinetic energy, Equation (8), calls for T ∼ [ d d τ x ( τ ) ] 2 , i.e.,

169 symbolic expressions are to be called upon. Each term is a complicated ω-dependent, one of these terms is the squared of Equation (10)! Justifiably we suppressed writing it down. Display of the kinetic and potential energies vs. ω and their intersection is shown in

Applying the FindRoot utility of Mathematica at about the guesstimated abscissa gives the exact root value of Equation (9), ω = 0.7992207. Reference [

To further our investigation, we compare the expansion coefficients of our approach vs. [

Having the values of expansion coefficients and ω we plot the t-dependent position x(t) given by Equation (5). This is shown in

Duffing equation, Equation (1) as discussed has an exact solution [

n | 1 | 3 | 5 | 7 | 9 | 11 | 13 | 15 | 17 | 19 | 21 | 23 | 25 | 27 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

a[n] | 1. | −2.3482 | 1.3617 | −1.4977 | 1.1270 | −1.0472 | 0.8581 | −0.7593 | 0.6390 | −0.5568 | 0.4727 | −0.4099 | 0.3489 | −0.3021 |

solution may be expressed in terms of the built in Mathematica library JacobiCN function,

x ( t ) = A J a c o b i C N [ 1 + ϵ A 2 t , ( ϵ A 2 2 ( 1 + ϵ A 2 ) ) 2 ] .

in Green. Although our numeric computation applying Power Series includes “only” 14-terms its accuracy is as close as the exact solution. We further our investigation applying Mathematica NDSolve solving Equation (1) numerically. NDSolve uses one of the common numeric methods of Implicit Runge Kutta or StiffnessSwitching which is a combination of the explicit and implicit methods solving the equation. The numeric solution perfectly matches the exact solution. The Green curve is the exact solution, it indistinguishably overlaps with the numeric solution. The black curve is the Power Series solution.

Despite the fact that the number of terms in Equation (5) is limited to fourteen the agreement between the exact/numeric and approximate results are quite acceptable. As shown in

As shown in

In this subsection we investigate the validity of the Power Series Method for a modified Duffing equation by replacing the ϵx^{3}(t) with ϵx^{4}(t). Being practical for the case on hand we consider { x [ 0 ] , ϵ } = { 0.5 , 0.5 } . By trial and error, we set the number of the terms, N, in Equation (5) to 20. Based on what we have discussed the non-zero terms limits the contributing terms to eleven. For the quartic term

the maximum potential energy is, V max ( t = 0 ) = 1 2 A 2 + 1 5 ϵ A 5 . Plotting the kinetic

and potential energies vs. ω their intersection abscissa according to Equation (9) yields ω. This is shown in

Numeric values of

Plot of Equation (5) is shown in

n | 1 | 3 | 5 | 7 | 9 | 11 | 13 | 15 | 17 | 19 |
---|---|---|---|---|---|---|---|---|---|---|

a[n] | 0.5 | −1.0203 | 0.0681 | −0.0745 | 0.0384 | −0.0119 | 0.0047 | −0.0037 | 0.0001 | −0.0012 |

In this subsection similar to Sect. 3.2 we investigate the validity of the Power Series Method for solving a modified Duffing equation with ϵx^{5}(t). We consider, { A , ϵ } = { 0.5 , 0.5 } . With trial and error, we set the number of the terms, N, in Equation (5) to 20. The maximum potential energy for the case on hand is

V max ( t = 0 ) = 1 2 A 2 + 1 6 ϵ A 6 . We skip displaying the kinetic and potential energies

vs. ω as is similar to

The numeric solution is the output of NDSolve. The equation on hand likewise to the previous case, Sect. 3.2, has no exact solution.

We applied a Computer Algebra System (CAS), Mathematica examining its applicability to solving Duffing equation utilizing the Power Series Method. We have shown with Mathematica we are able producing the numeric results of [

Moreover, we extend the calculation to symbolic domain; this is not included in [

For Duffing equation computation power of a laptop with a double processor prevents increasing the number of the terms evaluating Equation (5) to achieve accuracies beyond what is reported. Equation (7) is a recursive relationship between the expansion coefficients embodying b[n]. In general we write

{ ∑ n = 1 N a [ n ] τ n − 1 } m = ∑ n = 1 N b [ n ] τ n − 1 , for m = 2 , 3 , 4 , ⋯ ; for Duffing case, m = 3. We

notice replacing N = 28 with 30 drastically increases the number of b[n] greatly elongating the computation run-time crashing the computer. Number of terms in [

Extension of the current investigation. A curious reader may inquire: “What are the implications to numeric and symbolic solutions of the Duffing, modified Duffing equations embodying quartic and quintic terms if the second initial condition is replaced with, x ˙ [ t = 0 ] ≠ 0 ?” Quantification of the answer would require accessibility to a much more powerful computer as it would encounter lengthy symbolic calculation.

A note to the readers: Figures and calculations in this article are carried out utilizing Mathematica resources [

Motivation of this investigation stems from the questions that the author put to the presenter of “Power Series solution for strongly non-linear conservative oscillators” at “6^{th} International Physics Conference”, Athens/Greece, July 2018. One of the questions was “...has the Power Series Method been tested for nonlinear terms higher than three?” the answer was negative. The author initiated the investigation crafting this note. The author thanks Professor Mazen I. Qaisi for private communications and appreciates the support staff of Wolfram Research Inc. for consulting the issues concerning the numeric algorithm used in NDSolve.

The authors declare no conflicts of interest regarding the publication of this paper.

Sarafian, H. (2018) Sinusoidal Time-Dependent Power Series Solution to Modified Duffing Equation. American Journal of Computational Mathematics, 8, 259-268. https://doi.org/10.4236/ajcm.2018.83021