A simple cubication method for approximate solution of nonlinear Hamiltonian oscillators

Abstract

A new cubication method is proposed for periodic solution of nonlinear Hamiltonian oscillators. The method is formulated based on quasi-static equilibrium of the original oscillator and the undamped cubic Duffing oscillator. The cubication constants derived from the present cubication method are always based on elementary functions and are simpler than the constants derived by other cubication methods. The present method was verified using three common examples of strongly nonlinear oscillators and was found to give reasonably accurate results. The method can be used to introduce nonlinear oscillators in relevant undergraduate physics and mechanics courses.

Keywords

Hamiltonian system Duffing oscillator quasi-static equilibrium cubication Jacobi elliptic function elliptic integral

Introduction

Hamiltonian or conservative systems are freely oscillating systems in which the total energy is conserved and there is no dissipation of energy. Hamiltonian systems are typical of many mechanical systems such as vibration of the simple pendulum,^1,2 mass-spring system,^2,3 slider-crank mechanism,² shafts,³ beams,³ marine risers/drill string^4,5 and micro-electromechanical devices.⁶ The vibration response of practical Hamiltonian systems is usually nonlinear and the governing differential equations are generally difficult to solve. This makes their dynamic analysis quite challenging. As a result, the nonlinear vibrations of Hamiltonian systems are not usually included in undergraduate physics and mechanics courses. Even when included, the solutions provided for the nonlinear vibration models are based on classical perturbation methods (e.g. see the classic undergraduate mechanical vibration text by Rao³) that are obviously limited to small-amplitude vibrations and weakly nonlinear oscillators. To address this challenge, simple and reasonably accurate solution schemes that can be understood and applied by undergraduate students without much algebraic manipulation or complex mathematical techniques are required.

An interesting Hamiltonian system that is well recognized amongst nonlinear dynamicist is the undamped cubic Duffing oscillator. Its motion is described mathematically as

\ddot{u} + α_{1} u + α_{3} u^{3} = 0

(1)

with initial conditions

u (0) = A

and

\dot{u} (0) = 0

, where

f (u) = α_{1} u + α_{3} u^{3}

is the nonlinear restoring force, and

α_{1}

and

α_{3}

are the linear and nonlinear stiffness constants, respectively.

It has been long established that equation (1) has a closed-form solution in terms of the Jacobi elliptic cn function as shown¹

u (t) = A c n (ψ t, m)

(2)

where

ψ = \sqrt{α_{1} + α_{3} A^{2}}

m = \frac{α_{3} A^{2}}{2 (α_{1} + α_{3} A^{2})}

is the modulus of the elliptic function, and

K (m) = \int_{0}^{π / 2} \frac{d θ}{\sqrt{1 - m \sin^{2} θ}}

is the complete elliptic integral of the first kind. Hence, the time period for equation (1) is given as

T = 2 π / ω = 4 K (m) / ψ

(3)

The

K (m)

and

c n

functions are available in Mathematica and MatLab, but their series approximation is given in Appendices 1 and 2. Hence, equation (2) can be easily evaluated. However, for fundamentals on Jacobi elliptic functions, elliptic integrals and their application in the solution of nonlinear oscillators, the interested reader is encouraged to read Kovacic et al.⁷

The idea behind “cubication” methods is to represent the motion of a nonlinear oscillator with an equivalent motion in the form of equation (1) so that equations (2) and (3) will become approximate solutions to the original nonlinear oscillator. Hence, a number of cubication methods have been proposed^1,8–10 to provide approximate solutions to nonlinear Hamiltonian systems. One way to determine the equivalent cubication oscillator is to get the Taylor series for the nonlinear restoring force $f (u)$ and truncate at the cubic term. Since Taylor series is an infinite series, truncating at the cubic term is accurate for only small amplitudes (i.e. $A < 1$ ) and the error in this approximation grows unbounded as the amplitude increases. Belendez et al.¹ reasoned that the weakness of the Taylor series method of cubication lies in the fact the constants $α_{1}$ and $α_{3}$ are independent of the oscillation amplitude, $A$ . Hence, they proposed a cubication method based on Chebyshev polynomial in which $α_{1}$ and $α_{3}$ are functions of $A$ . The Chebyshev polynomial was also applied in Elıas-Zuniga and Martınez-Romero¹⁰ to formulate an enhanced cubication method. The main advantage of the Chebyshev series expansion is that it converges more quickly than Taylor series¹¹ and it produces a series approximation in cases where Taylor series cannot produce an approximation.¹ The Chebyshev cubication methods in literature^1,10 produced approximate expressions for the frequency–amplitude response that were found to be in good agreement with exact solutions for some examples of oscillators with odd nonlinearity. We note that cubication methods are suited for odd nonlinear Hamiltonian oscillators because the restoring force of equation (1) is an odd function.

The challenge with the Chebyshev cubication methods is that deriving $α_{1} (A)$ and $α_{3} (A)$ for some oscillators is too complicated for undergraduate courses. For example, in the case of the nonlinear simple pendulum where the restoring force is $f (u) = \sin u$ , $α_{1} (A)$ and $α_{3} (A)$ were derived in terms of special functions, i.e. first- and third-order Bessel functions of the first kind.¹ In the present study, a simple cubication method that is based on quasi-static equilibrium conditions is proposed to formulate approximate periodic solutions to odd nonlinear Hamiltonian oscillators. The expressions for $α_{1} (A)$ and $α_{3} (A)$ derived based on the present method are simple and always in terms of elementary functions. The present cubication method was verified using three examples of strong nonlinear oscillators and found to give good accuracy compared to standard numerical solutions. The purpose of this paper is to present a solution method that can be used to introduce nonlinear Hamiltonian oscillators in relevant undergraduate physics and mechanics courses.

Present cubication method

Let us consider a nonlinear oscillator with odd restoring force ( $g (\pm u) = \pm g (u)$ ) given as

\ddot{u} + g (u) = 0

(4)

where

u (0) = A

and

\dot{u} (0) = 0

. The present method is based on the idea that the original (equation (4)) and equivalent (equation (1)) oscillators must have the same restoring force and potential energy at the amplitude. These quasi-static equilibrium conditions require that

for force equilibrium condition, $f (A) = g (A)$ and

for potential energy equilibrium condition, $\int_{0}^{A} f (u) d u = \int_{0}^{A} g (u) d u$ .

The above assumptions are reasonable from a physical viewpoint because the restoring force of the equivalent cubic oscillator is meant to approximate that of the original nonlinear oscillator. Therefore, the restoring force and potential energy of both oscillators should be similar. The above quasi-static equilibrium conditions assume that both oscillators have identical restoring force and potential energy at the critical point where the displacement is maximum. If these quasi-static assumptions work well then the restoring force of the equivalent and original oscillator would be comparable. Evidence of such agreement in the restoring force prediction of both oscillators is shown in the next section.

From equation (1), we get $f (A) = α_{1} A + α_{3} A^{3}$ and $\int_{0}^{A} f (u) d u = \frac{1}{2} α_{1} A^{2} + \frac{1}{4} α_{3} A^{4}$ . Therefore, substituting these results into the above conditions we have

α_{1} A + α_{3} A^{3} = g (A)

(5a)

2 α_{1} A^{2} + α_{3} A^{4} = 4 \int_{0}^{A} g (u) d u

(5b)

Solving equations (5a) and (5b) simultaneously gives

α_{1} (A) = \frac{4 \int_{0}^{A} g (u) d u - A g (A)}{A^{2}}

(6a)

α_{3} (A) = \frac{2 A g (A) - 4 \int_{0}^{A} g (u) d u}{A^{4}}

(6b)

Equations (6a) and (6b) give the cubication constants in terms of $A$ as in the case of the Chebyshev cubication methods. However, the present cubication method is much simpler than the Chebyshev cubication methods.

Periodic solutions for some oscillators with strong nonlinearity

Three examples of strong nonlinear Hamiltonian vibrations of mechanical systems were used to confirm the accuracy of the present cubication method. Equations (6a) and (6b) were applied to derive the cubication constants, while equations (2) and (3) were used to determine the periodic solutions, i.e. time period and oscillation history.

Cubic-quintic Duffing oscillator

The motion of the cubic-quintic oscillator can be expressed as¹²

\ddot{u} + k_{1} u + k_{3} u^{3} + k_{5} u^{5} = 0

(7)

where

g (u) = k_{1} u + k_{3} u^{3} + k_{5} u^{5}

. Therefore,

α_{1} = k_{1} - \frac{1}{3} k_{5} A^{4}

and

α_{3} = k_{3} + \frac{4}{3} k_{5} A^{2}

. A comparison of the periodic solution for the present cubication method and numerical integration results when

k_{1} = k_{3} = k_{5} = 1.0

and

A = 1.0

is shown in Figure 1(a) and the restoring force is shown in Figure 1(b). A variation of the time period with amplitude is presented in Table 1 and the corresponding error analysis is plotted in Figure 2. The percentage relative error converges to a maximum value of 3.54% as

A \to \infty

. These results show a good agreement between the present method and very accurate solutions obtained by numerical integration of equation (7) using the StiffnessSwitching method⁶ of the NDSolve function in Mathematica.

Figure 1.

(a) Periodic solution and (b) restoring force: cubic-quintic oscillator with $k_{1} = k_{3} = k_{5} = 1.0$ and $A = 1.0$ .

Figure 2.

Error analysis for the time period estimates of the present cubication method: cubic-quintic oscillator with $k_{1} = k_{3} = k_{5} = 1.0$ .

Table 1.

Comparison of time period estimates: cubic-quintic oscillator.

Amplitude	Present method	Numerical
0.01	6.28295	6.28284
0.1	6.25956	6.26044
0.5	5.67217	5.67793
1	4.08739	4.12391
5	0.316571	0.32756
10	0.080662	0.0835769
100	0.000811583	0.000841217
1000	8.11583 × 10^–6	8.41293 × 10^–6
10,000	8.11584 × 10^–8	8.41294 × 10^–8
100,000	8.11584 × 10^–10	8.41365 × 10^–10
1,000,000	8.11584 × 10^–12	8.41276 × 10^–12

Simple pendulum

The model for the nonlinear oscillations of an undamped simple pendulum can be written as^1,2

\ddot{u} + Ω \sin u = 0

(8)

where

g (u) = Ω \sin u

;

Ω = g / l

and

0 < A < π

. Therefore,

α_{1} = Ω [4 (1 - \cos A) - A \sin A] / A^{2}

and

α_{3} = Ω [2 A \sin A - 4 (1 - \cos A)] / A^{4}

. These expressions for

α_{1}

and

α_{3}

are simpler compared to the corresponding expressions in literature^1,10 which were derived in terms of Bessel functions. The periodic solutions for

A = π / 2

(i.e.

90 °

) and

Ω = 1.0

obtained using the present cubication method and numerical integration of equation (8) are shown in Figure 3(a), while the restoring force is shown in Figure 3(b). The periodic solution of the present method is indistinguishable from the numerical results. Also, the time period estimates for

0 < A < 180 °

when

Ω = 1.0

are given in Table 2, while the corresponding error estimates are plotted in Figure 4. It was observed that the maximum error in the time period estimate of the present method is less than 1% for

A \leq 135 °

, less than 5% for

A \leq 171.50 °

and less than 10% for

A \leq 179 °

. This is considered a good approximation in view of the simplicity of the present method.

Figure 3.

(a) Periodic solution and (b) restoring force: simple pendulum with $Ω = 1.0$ and $A = π / 2$ radians.

Figure 4.

Error analysis for the time period estimates of the present cubication method: simple pendulum with $Ω = 1.0$ and $0 < A \leq 150 °$ .

Table 2.

Comparison of time period estimates: simple pendulum.

Amplitude	Present method	Numerical
5°	6.28618	6.28882
10°	6.29517	6.29532
20°	6.33135	6.3305
30°	6.39248	6.39297
40°	6.47983	6.47999
50°	6.59527	6.59588
60°	6.74146	6.74251
70°	6.92194	6.9249
80°	7.1415	7.1473
90°	7.40656	7.41771
100°	7.72591	7.74211
110°	8.11183	8.14062
120°	8.58214	8.62608
130°	9.164	9.23684
135°	9.50961	9.60026
140°	9.90196	10.0218
150°	10.8775	11.072
160°	12.269	12.6137
170°	14.6304	15.3268
179°	22.1364	24.5133

Mass attached to two stretched elastic springs

The nonlinear vibration of a mass attached to two stretched elastic springs is governed by¹

\ddot{u} + ω_{0}^{2} (1 - \frac{λ}{\sqrt{1 + u^{2}}}) u = 0

(9)

where

g (u) = ω_{0}^{2} (1 - \frac{λ}{\sqrt{1 + u^{2}}}) u

ω_{0}

is a constant and

0 < λ \leq 1

. Therefore,

α_{1} = ω_{0}^{2} [1 + \frac{4 λ (1 - \sqrt{1 + A^{2}})}{A^{2}} + \frac{λ}{\sqrt{1 + A^{2}}}]

and

α_{3} = - \frac{2 λ ω_{0}^{2}}{A^{2}} [\frac{1}{\sqrt{1 + A^{2}}} + \frac{2 (1 - \sqrt{1 + A^{2}})}{A^{2}}]

. Again, these expressions for

α_{1}

and

α_{3}

are simpler compared to the corresponding expressions in Belendez et al.,¹ which were derived in terms of elliptic functions of the first and second kind. Figure 5(a) shows the periodic solution for

A = 1.0

and

λ = 1.0

obtained using the present cubication method and numerical integration of equation (9). Also, the restoring force is shown Figure 5(b). Since the response of this system depends on the parameter

λ

, the time period was investigated for different values of

λ

as shown in Table 3. The corresponding error analysis is plotted in Figure 6, which shows that the error in the present method is approximately zero when

A \to 0

and when

A \to \infty

. This observation agrees with the results in Belendez et al.¹

Figure 5.

(a) Periodic solution and (b) restoring force: mass attached to two stretched elastic springs with $λ = 1.0$ and $A = 1.0$ .

Figure 6.

Error analysis for the time period estimates of the present cubication method: mass attached to two stretched elastic springs.

Table 3.

Comparison of time period estimates: mass attached to two stretched elastic spring.

Amplitude	$λ = 0.25$		$λ = 0.75$		$λ = 1.00$
Amplitude	Present method	Numerical	Present method	Numerical	Present method	Numerical
0.01	7.25515	7.25491	12.5657	12.5629	1048.86	1049.52
0.1	7.2507	7.25331	12.4968	12.5027	105.215	105.192
0.5	7.15969	7.16047	11.2837	11.2707	22.5429	22.4447
1	6.99016	6.9854	9.66828	9.6254	13.2115	13.0676
2	6.7493	6.73423	8.12537	8.04819	9.24928	9.09865
5	6.50005	6.48262	7.01454	6.94942	7.32502	7.22611
10	6.39646	6.38409	6.64368	6.60211	6.77918	6.72036
100	6.29492	6.29328	6.31861	6.31347	6.33055	6.32357
1000	6.28436	6.28419	6.28672	6.28622	6.2879	6.28721
10,000	6.28354	6.28324	6.28354	6.28356	6.28366	6.28362
100,000	6.2832	6.28775	6.28322	6.28339	6.28323	6.28334
1,000,000	6.28319	6.28315	6.28319	6.28319	6.28319	6.28338

Note: As shown in Figure 6, the maximum relative errors of the present method for $λ = 0.25$ , $0.75$ and 1.00 are 0.269%, 0.959% and 1.656%, respectively. The maximum errors were found to occur in the amplitude range of $1 < A < 10$ .

Concluding remarks

This paper presents a new cubication method that is based on quasi-static equilibrium considerations. The cubication constants derived are simple and straightforward compared to existing cubication methods based on Chebyshev series,^1,10 weighted mean-square⁸ and harmonic balance principles.⁹ Furthermore, the cubication constants of the present method are always expressed as elementary functions provided the restoring force is based on elementary functions.

The present method can be used to introduce nonlinear vibrations of Hamiltonian oscillators in undergraduate physics and mechanics courses due to its simplicity and accuracy for a wide range of problems. It can also add pedagogical value to undergrads by introducing the Jacobi elliptic function and the elliptic integral of the first kind in an appealing manner. This is because the present method applies these special functions in a simple scheme that does not make them appear frightening.

An important quality of the present cubication method is that it is formulated based on the mechanics (i.e. restoring force and potential energy) of the original nonlinear oscillator and it is not just a purely mathematical exercise. Hence, the present method can add to understanding of how the mechanics of nonlinear oscillators can be applied in obtaining their periodic solutions. This quality is also shared by some energy-based methods for periodic solution of nonlinear oscillators, e.g. the Energy Balance Method.¹³ Furthermore, the present cubication method can be used for qualitative analysis of the original nonlinear oscillator.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Belendez

Alvarez

Fernandez

, et al. Cubication of conservative nonlinear oscillators. Eur J Phys 2009; 30: 973–981.

Fidlin

Nonlinear oscillations in mechanical engineering. New York: Springer, 2005.

Rao

SS.

Mechanical vibrations. 5th ed. Upper Saddle River: Pearson Education Inc., 2011.

Major

Big-Alabo

Odi-Owei

Effect of drill string rotation on the dynamic response of drilling risers. Int J Appl Mech Eng 2015; 20: 503–516.

Kaewunruen

Chiravatchradej

Chucheepsakul

Nonlinear free vibration of marine risers/pipes transporting fluids. Ocean Eng 2005; 32: 417–440.

Big-Alabo

Continuous piecewise linearization method for approximate periodic solution of the relativistic oscillator.

Int J Mech Eng Educ. 2020; 48: 178–194.

Kovacic

Cveticanin

Zukovic

, et al. Jacobi elliptic functions: a review of nonlinear oscillatory application problems. J Sound Vib 2016; 380: 1–36.

Bravo Yutse

Martin Sanchez

A weighted mean-square method of “cubication” for nonlinear oscillators. J Sound Vib 1989; 134: 423–433.

Bravo Yutse

“Cubication” for nonlinear oscillators using the principle of harmonic balance. Int J Non-Linear Mech 1990; 27: 347–356.

10.

Elıas-Zuniga

Martınez-Romero

Accurate solutions of conservative nonlinear oscillators by the enhanced cubication method.

Math Probl Eng 2013: 1–9. Article 842423. DOI:10.1155/2013/842423

11.

Elias-Zuniga

“Quintication” method to obtain approximate analytical solutions of non-linear oscillators. Appl Math Comput 2014; 243: 849–855.

12.

Big-Alabo

Periodic solutions of Duffing-type oscillators using continuous piecewise linearization method. Mech Eng Res 2018; 8: 41–52.

13.

JH.

Preliminary report on the energy balance for nonlinear oscillations. Mech Res Commun 2002; 29: 107–111.