Locally exact discretization method for nonlinear oscillation systems

Abstract

Nonlinear oscillations are ubiquitous in various fields of engineering and applied sciences. In engineering, vibration isolators such as absorbers and dampers have been used to reduce damages of most undesired vibrations to a structure or a system. One important class of dampers is the pendulum tuned mass damper. Since it is difficult to solve nonlinear oscillation equations analytically, it is highly demanded to design numerical methods that preserve as much as possible intrinsic properties of the original system. To this end, we propose a locally exact modified midpoint rule with a parameter ω by combining ideas of locally exact discretization and symplecticness. We prove that it is convergent of order two for any ω as its classical partner; however, its local error accuracy can be remarkably improved with a proper choice of ω. A numerical example is provided to highlight the accuracy and effectiveness of the new method. Finally, effective implementation of the new method is also discussed based on approximate analytical solution.

Keywords

Tuned mass damper exact solution Hamiltonian systems modified midpoint rule simple pendulum

Introduction

Nonlinear oscillations are ubiquitous in various fields of engineering and applied sciences,¹ for example swinging pendulum, the beating of the human heart, business cycles in economics, etc. Study and simulation of nonlinear oscillation systems are important. In engineering, to reduce damages of most undesired vibrations to a structure or a system, vibration isolators such as absorbers and dampers have been used widely, for example the designer of the Taipei 101 building designed a 660-tonne steel pendulum that serves as a tuned mass damper (TMD) to withstand the typhoon winds and earthquake tremors that are common in the area east of Taiwan. In order to understand the behavior of a pendulum TMD in detail in practice, it is often needed to solve efficiently a mathematical modelling of nonlinear pendulum oscillation. Due to Newton’s second law, most motion of undamped oscillations in scientific and engineering applications are modeled by the second differential equation

\ddot{q} = f (q)

(1)

where

f : ℝ^{n} \to ℝ^{n}

is assumed to be sufficiently smooth so that the equation (1) has a unique solution under appropriate initial conditions.

It is well-known that it is difficult to obtain explicit solution for nonlinear equation (1) analytically, therefore, to find accurate enough numerical solution is desirable at proper cost with small step size. However, the dynamics of conventional numerical methods such as Runge-Kutta method, multi-step method, predictor-corrector method and others may not be as rich as that of its continuous counterparts.^2,3 This type of inaccuracy cannot be improved by reducing the step size of numerical methods. Therefore, numerical methods that preserve as much as possible intrinsic properties of the original dynamical system for all step sizes or relatively big ones are highly demanded. Exponential method which exactly integrates the linearized equations for semilinear or nonlinear equations has received much attention over last few decades.⁴ In this direction, some instances of this type of integrators are methods known as exponential fitting,^5–7 exponential integrators,^8,9 local linearization^10,11 and locally exact modification of conventional numerical methods.^12,13

By introducing $p = \dot{q},$ equation (1) could be written in separable Hamiltonian systems with $f (q) = \nabla F (q)$ . Two important features of Hamiltonian systems include that the flow is a symplectic transform and the Hamiltonian function (sometimes called energy) is an invariant. Numerical method that preserves symplectic structure of Hamiltonian systems is called symplectic, and numerical method which preserves Hamiltonian function or energy is called energy-preserving. Symplectic method, energy-preserving method and other structure-preserving methods have been wildly recognized (see review papers^14–16 and the references therein) that these methods lead to an improved long-time behavior over general traditional methods.

Recently, Cieśliński and Ratkiewicz¹³ successfully present locally exact discrete gradient for equation (1), where locally exact modification means that the linearization of modified discrete gradient scheme around some point gives exact discretization of the linearization of equation (1) around the same point. Analysis and numerical experiments show that locally exact discrete gradient method has attractive advantage while retaining energy conserved. On the other hand, it is well-known that there exists no numerical method that can be simultaneously symplectic and energy-preserving, otherwise it is a time-transformation of the exact flow of Hamiltonian systems.^17,18 Therefore, construction of symplectic algorithm is of great interest in numerical integration of Hamiltonian systems. In addition, by a backward error analysis, symplectic algorithm could exhibit bounded energy error for exponentially long times (see, e.g. Hairer et al.¹⁹).

Motivated by the importance of symplecticness for Hamiltonian systems and the advantage of locally exact discretization suggested by Cieśliński and Ratkiewicz,¹³ we hope to discuss to what extent and how locally exact discretization compete with symplectic scheme for the considered equation (1), and how locally exact modified symplectic scheme provides improved numerical simulation for nonlinear oscillation systems. To this end, our work here is to present locally exact symplectic algorithm. We take the implicit midpoint rule as an example and construct its locally exact modification while keeping symplectic condition simultaneously, aiming to provide better numerical simulation for nonlinear oscillation systems. Before we close this introduction, it is worth mentioning that apart from numerical methods for solving equation (1), recent semi-analytical methods that focus on frequency-amplitude estimation and active vibration control^20–24 are interesting, which give approximate analytical solution of equation (1).

This paper is organized as follows. In section Preliminary, basic knowledge about numerical integration of Hamiltonian systems are briefly reviewed. In section Locally exact modified midpoint rule, the process to construct locally exact modified midpoint rule is proposed and convergence of this method is analyzed. In section Numerical Example, numerical experiment is presented to show that the proposed method possesses expected accuracy and better numerical simulation for simple pendulum equation. Finally, this paper ends with some concluding remarks and discussions.

Preliminary

For completeness, we firstly present some basic knowledge about Hamiltonian systems and their relevant numerical discretization methods, which could be found in several classical monographs, for example see literature.^19,25–27

A d-dimensional autonomous Hamiltonian system consists of 2d first-order differential equations

\dot{y} = F (y) = J^{- 1} \nabla H (y), J = [\begin{matrix} 0 & - I_{d} \\ I_{d} & 0 \end{matrix}]

(2)

where I_d is d × d identity matrix. Denote

y = {(q_{1}, \dots, q_{d}, p_{1}, \dots, p_{d})}^{T}

, equation (2) could be written in component form

\dot{q_{k}} = \frac{\partial H}{\partial p_{k}}, \dot{p_{k}} = - \frac{\partial H}{\partial q_{k}}, k = 1, \dots, d

(3)

where H(q, p) is called Hamiltonian function. When

H (q, p) = T (q) + V (p)

it is called separable Hamiltonian systems.

The flow $φ_{t}$ of the Hamiltonian systems (2) is a symplectic transform, i.e. the Jacobian derivative of $φ_{t}$ with respect to y satisfies the matrix identity

φ_{t}' {(y)}^{T} J φ_{t}' (y) = J

(4)

A desirable property of a numerical method defined by $y_{n + 1} = ψ_{h} (y_{n})$ for equation (2) is that in addition to providing an accurate enough approximation of the exact solution for small enough step size h, it preserves the qualitative property of the continuous flow $φ_{t}$ such as symplecticness. This motivates the following definition by mimicking the equality (4).

Definition 2.1. A numerical method $y_{n + 1} = ψ_{h} (y_{n})$ is symplectic for Hamiltonian systems equation (2) if the Jacobian matrix of the numerical flow satisfies

ψ_{h}' {(y_{n})}^{T} J ψ_{h}' (y_{n}) = J

(5)

The well-known approaches to construct symplectic method include generating function theory,²⁸ variational method,^29,30 the class of Runge-Kutta (RK) method (including partitioned Runge-Kutta method and Runge-Kutta-Nyström method), splitting and composition method. In this paper, we mainly focus on RK method.

An s-stage RK method for solving (2) with the initial condition $y (t_{0}) = y_{0}$ is a one-step method defined by

\begin{array}{l} Y_{i} = y_{0} + h \sum_{j = 1}^{s} a_{i j} F (Y_{j}), i = 1, \dots, s \\ y_{1} = y_{0} + h \sum_{i = 1}^{s} b_{i} F (Y_{i}) \end{array}

(6)

where h is the step size. Usually the coefficients of the above-defined RK method are displayed by the following Butcher tabular

\begin{array}{l} c_{1} & a_{11} & \dots & a_{l s} \\ ⋮ & ⋮ & ⋮ \\ c_{s} & a_{s 1} & \dots & a_{s s} \\ b_{1} & \dots & b_{s} \end{array}

A characteristic for an RK method is symplectic if its coefficients satisfy

b_{i} a_{i j} + b_{j} a_{j i} = b_{i} b_{j}, i, j = 1, \dots, s

(7)

Classical midpoint rule is a one-stage RK method with Butcher tabular

\begin{matrix} \frac{1}{2} & \frac{1}{2} \\ 1 \end{matrix}

(8)

It is easy to verify that the coefficients in (8) satisfy that of equation (7), therefore, the midpoint rule is a symplectic method.

In fact, if we adjust the coefficients in (8) by

\begin{matrix} c_{1} & \frac{a}{2} \\ a \end{matrix}

(9)

then it is easy to obtain that the above RK method also satisfies (7) for any meaningful parameter a and consistent c₁.

In the next section we will determine the parameter a in (9) by forcing this method solves exactly the linearized equation of (1). We call it locally exact modified midpoint rule.

Locally exact modified midpoint rule

Now we are in the position to present our locally exact modified midpoint rule in the following steps.

Step 1 By taking linear Taylor expansion of f(q) around some point $\bar{q},$ we get the linearized equation (1)

\ddot{q} = f' (\bar{q}) q + f (\bar{q}) - f' (\bar{q}) \bar{q}

(10)

By introducing $p = \dot{q}, w^{2} \equiv f' (\bar{q}), g \equiv f (\bar{q}) - f' (\bar{q}) \bar{q}$ , then equation (10) could be written by linear Hamiltonian equations

\begin{array}{l} \dot{q} = p, \\ \dot{p} = f' (\bar{q}) q + g \equiv w^{2} q + g \end{array}

(11)

The exact solutions of equation (11) with initial conditions $q (0) = q_{0}, p (0) = p_{0}$ are

\begin{array}{l} q (t) = \frac{1}{2} (e^{w t} + e^{- w t}) q_{0} + \frac{1}{2 w} (e^{w t} - e^{- w t}) p_{0} + \frac{g}{2 w^{2}} (e^{w t} + e^{- w t} - 2), \\ p (t) = \frac{w}{2} (e^{w t} - e^{- w t}) q_{0} + \frac{1}{2} (e^{w t} + e^{- w t}) p_{0} + \frac{g}{2 w} (e^{w t} - e^{- w t}) \end{array}

(12)

for

w \neq 0.

w = 0,

we take the limit as exact solutions.

Step 2 By applying the modified midpoint rule (9) with the parameter a to equation (11), we get the following numerical solutions for one step

\begin{array}{l} q_{1} = q_{0} + h a \frac{p_{1} + p_{0}}{2}, \\ p_{1} = p_{0} + h a (w^{2} \frac{q_{1} + q_{0}}{2} + g) \end{array}

(13)

Therefore, we get

\begin{array}{l} q_{1} = \frac{4 + a^{2} w^{2} h^{2}}{4 - a^{2} w^{2} h^{2}} q_{0} + \frac{4 h a}{4 - a^{2} w^{2} h^{2}} p_{0} + \frac{2 a^{2} g h^{2}}{4 - a^{2} w^{2} h^{2}}, \\ p_{1} = \frac{4 a w^{2} h}{4 - a^{2} w^{2} h^{2}} q_{0} + \frac{4 + a^{2} w^{2} h^{2}}{4 - a^{2} w^{2} h^{2}} p_{0} + \frac{4 h a g}{4 - a^{2} w^{2} h^{2}} \end{array}

(14)

Step 3 By making the numerical solutions q₁, p₁ obtained by equation (14) equal to the exact solutions of (12) at $t = h,$ we derive the unknown parameter a in modified midpoint rule (9), which is

\begin{array}{l} a = \frac{2}{w h} \tanh (\frac{w h}{2}) \end{array}

(15)

where

w^{2} = f' (\bar{q}) .

When

w h \to 0,

we take

a = 1.

Step 4 Once we determine the parameter a in Butcher tabula (9), we get the modified midpoint rule for the original nonlinear equation (1)

\begin{array}{l} q_{n + 1} = q_{n} + \frac{2}{w} \tanh (\frac{w h}{2}) \frac{p_{n + 1} + p_{n}}{2}, \\ p_{n + 1} = p_{n} + \frac{2}{w} \tanh (\frac{w h}{2}) f (\frac{q_{n + 1} + q_{n}}{2}) \end{array}

(16)

where

w^{2} = f' (\bar{q})

, and

\bar{q}

is some point which could be decided by different strategies. We will discuss this point in the last section of the paper.

Remark 3.1. The modified midpoint rule (16) is consistent with the locally exact implicit midpoint rule in Cieśliński¹² (see (3.18) therein). Cieśliński¹² presents it as locally exact implicit midpoint rule by making its linearization around some point is identical to the exact discretization of the equation (1) linearized at the same point. However, we here derive (16) in a slightly different approach that we determine the parameter a in (9) by forcing this method solves exactly the linearized equation of (1).

Theorem 3.1. For smooth f(q) in equation (1) with bounded derivatives up to the third order, the numerical scheme (16) has the second-order convergence. Particularly, when $w^{2} = f' (q_{n})$ , the local error accuracy for q is the third order.

Proof: We prove the conclusion for one-step approximation from (q₀, p₀) to $(q_{1}, p_{1}) .$ Firstly, we expand exact solutions $q (t_{0} + h), p (t_{0} + h)$ at (q₀, p₀) in Taylor series

q (t_{0} + h) = q_{0} + h p_{0} + \frac{h^{2}}{2} f (q_{0}) + \frac{h^{3}}{6} f' (q_{0}) p_{0} + O (h^{4})

(17)

p (t_{0} + h) = p_{0} + h f (q_{0}) + \frac{h^{2}}{2} f' (q_{0}) p_{0} + \frac{h^{3}}{6} (f' (q_{0}) f (q_{0}) + f ″ (q_{0}) p_{0}^{2}) + O (h^{4})

(18)

Secondly, we expand numerical solutions q₁, p₁ at (q₀, p₀) in the following

q_{1} = q_{0} + a p_{0} + \frac{1}{2} a^{2} f (q_{0}) + \frac{1}{4} f' (q_{0}) p_{0} a^{3} + O (a^{4})

(19)

p_{1} = p_{0} + a f (q_{0}) + \frac{1}{2} a^{2} f' (q_{0}) p_{0} + (\frac{1}{4} f' (q_{0}) f (q_{0}) + \frac{1}{8} f ″ (q_{0}) p_{0}^{2}) a^{3} + O (a^{4})

(20)

Using the approximation $a = h - \frac{1}{12} w^{2} h^{3} + O (h^{5})$ equations (19) and (20) could be expressed by

q_{1} = q_{0} + h p_{0} + \frac{1}{2} h^{2} f (q_{0}) + (\frac{1}{4} f' (q_{0}) p_{0} - \frac{1}{12} w^{2} p_{0}) h^{3} + O (h^{4})

(21)

\begin{array}{l} p_{1} = p_{0} + h f (q_{0}) + \frac{1}{2} h^{2} f' (q_{0}) p_{0} + (\frac{1}{4} f' (q_{0}) f (q_{0}) + \frac{1}{8} f ″ (q_{0}) p_{0}^{2} - \frac{1}{12} w^{2} f (q_{0})) h^{3} + O (h^{4}) \end{array}

(22)

Subtracting equations (21) and (22) from equations (17) and (18), we arrive at

\begin{array}{l} q (t + h) - q_{1} = (\frac{1}{12} w^{2} p_{0} - \frac{1}{12} f' (q_{0}) p_{0}) h^{3} + O (h^{4}), \\ p (t + h) - p_{1} = (\frac{1}{12} w^{2} f (q_{0}) - \frac{1}{12} f' (q_{0}) f (q_{0}) + \frac{1}{24} f ″ (q_{0}) p_{0}^{2}) h^{3} + O (h^{4}) \end{array}

Therefore, it suggests that the accuracy of modified midpoint rule is the second order for any $ω .$ Specially, for the component q, the local error accuracy is the third order when $w^{2} = f' (q_{0})$ , which is higher than classical midpoint rule.□

Numerical example

In this section, we compare our modified midpoint rule with classical midpoint rule to show that the new modified midpoint rule possesses expected accuracy of the second order (see Table 1), but produces smaller global errors than classical midpoint rule (see Figure 1).

Figure 1.

Global errors of θ_n by modified midpoint rule (16) with red line and classical midpoint rule by blue line for $t \in [0, 1000], h = 0.1, ω_{0} = 1$ .

Table 1.

Global errors $E_{h}^{θ}$ between θ_n and $θ (t)$ at t = 10 and numerical orders.

	Midpoint rule (8)		Modified midpoint rule (16) ( $ω^{2} = f' (q_{n})$ )
h	$E_{h}^{θ} = \| θ_{n} - θ (10) \|$	Order	$E_{h}^{θ} = \| θ_{n} - θ (10) \|$	Order
$2^{- 3}$	6.041E-04		2.715E-07
$2^{- 4}$	1.525E-04	1.9865	0.689E-07	1.9786
$2^{- 5}$	0.382E-04	1.9966	0.173E-07	1.9919
$2^{- 6}$	0.096E-04	1.9992	0.043E-07	1.9966
$2^{- 7}$	0.024E-04	1.9998	0.011E-07	1.9986
$2^{- 8}$	0.006E-04	1.9999	0.003E-07	1.9999

To this end, we take one of the most popular nonlinear oscillation systems—simple pendulum—as an example, which is written as

l \frac{d^{2} θ}{d t^{2}} = - g \sin θ

(23)

where l is the length of thin and massless rigid rod, g is the acceleration of gravity,

θ (t)

denotes the angular displacement of the motion taking place in a vertical plane and t is the time. With initial conditions

θ (0) = θ_{0}, \dot{θ} (0) = 0,

equation (23) has exact analytical solution³¹

θ (t) = 2 \arcsin {\sin \frac{θ_{0}}{2} s n [K (\sin^{2} \frac{θ_{0}}{2}) - ω_{0} t; \sin^{2} \frac{θ_{0}}{2}]}

(24)

where

ω_{0} = \sqrt{\frac{g}{l}}, K (\cdot)

is complete elliptic integral of the first kind, and

s n [\cdot; \cdot]

is the Jacobi elliptic function.

Global errors $E_{h}^{θ} = | θ_{n} - θ (t_{n}) |$ between numerical solution and exact solution corresponding to classical midpoint rule and the modified midpoint rule with six small step sizes are shown in Table 1, where initial values are chosen as $θ_{0} = 5^{°}, \dot{θ} (0) = 0$ and $ω_{0} = 1$ . In addition, we define the order of a numerical method by

order \approx \frac{l n (E_{h}^{θ} / E_{\frac{h}{2}}^{θ})}{l n (2)}

From the third and fifth columns of Table 1, it is clearly seen that both classical midpoint rule and the modified midpoint rule possess the expected second order of convergence. Furthermore, the second and fourth columns of Table 1 show that modified midpoint rule obviously produces smaller errors than classical midpoint rule.

Next, Figure 1 shows global errors for $| θ_{n} - θ (t_{n}) |$ under different initial amplitude $θ_{0} .$ From every subplot, it is easy to see again that modified midpoint rule method obviously produces smaller errors than classical midpoint rule, even for large initial amplitude $θ_{0} .$

Conclusions and discussions

Exact discretization of linear differential equations is natural and have the very attractive property of leading to exact approximate solution. Symplectic algorithms are widely recognized for the numerical integration of Hamiltonian systems. In this paper, by combining these two effective ideas, we derive a locally exact modified midpoint rule (16) for the conservative nonlinear oscillation systems (1). We show that the modified midpoint rule (16) with any ω always has the second-order convergence. Particularly, if $w^{2} = f' (q_{n}),$ then the scheme (16) has improved local truncation error for q compared with classical midpoint rule. Numerical experiment verifies that the locally exact discretization method is extremely effective for nonlinear oscillation systems, and it provides a universal way to optimize a nonlinear oscillator including pendulum TMD system.

In order to implement the modified midpoint rule, one needs to determine a functional parameter $ω .$ Noting that $ω^{2} = f' (\bar{q}),$ some possible choices of $\bar{q}$ are initial values q₀, numerical approximation q_n at $t_{n} = t_{0} + n h,$ or an average numerical approximation $\frac{q_{n} + q_{n + 1}}{2} .$ In the numerical experiments of this paper, we choose $\bar{q} = q_{n} .$ However, a fixed or already known $\bar{q}$ is preferable in order to reduce computational costs on the number of nonlinear iteration. Therefore, it is very interesting to use approximate analytical solution of q as $\bar{q} .$ For the nonlinear oscillation systems (1) are conservative Hamiltonian systems, it is much convenient to obtain approximate analytical solution by He’s Hamiltonian approach.²⁰ For equation (1) with initial conditions $q (0) = q_{0}, \dot{q} (0) = 0,$ Hamiltonian approach gives an approximate analytical solution by $q (t) = q_{0} \cos (\tilde{ω} t),$ where $\tilde{ω}$ is determined by

\frac{\partial}{\partial q_{0}} (\frac{1}{\partial (\frac{1}{\tilde{ω}})} \int_{0}^{\frac{T}{4}} [\frac{1}{2} {\dot{q}}^{2} + F (q)] d t) = 0

with

\nabla F (q) = f (q), T = \frac{2 π /}{\tilde{ω}} .

Once we get approximate analytical solution

q (t),

it is easy to set

\bar{q} = q (t_{n})

in (16) while retaining the second-order convergence.

Last but not least, it is noted that there are some interesting issues related to this work which will be further considered.

In this paper, we consider efficient numerical computation of undamped nonlinear oscillation systems. It would be interesting to consider the damped case in the same spirit of this paper.

The modified midpoint rule is an implicit numerical method as the classical one. Since the problem we consider are separable Hamiltonian systems, it will be interesting to consider modification of explicit symplectic methods, such as Stömer-Verlet method.

Footnotes

Acknowledgements

The author thanks the anonymous referees and the Editor for their valuable comments and suggestions. Part of the work was done when the author was visiting the State Key Laboratory of Scientific and Engineering Computing, the Academy of Mathematics and Systems Science of the Chinese Academy of Sciences. The author thanks this institute for its support and hospitality.

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (No. 11761033).

References

Gitterman

The Noisy Oscillator: The First Hundred Years, From Einstein Until Now. Hackensack, NJ: World Scientific Publishing Co. Pte. Ltd, 2005.

Ascher

Reich

The midpoint scheme and variants for Hamiltonian systems: advantages and pitfalls. SIAM J Sci Comput 1999; 21: 1045–1065.

Stewart

Numerical methods: warning-handle with care!

Nature 1992; 355: 16.

Pope

An exponential method of numerical integration of ordinary differential equations. Commun Acm 1963; 6: 491–493.

Vanden Berghe

De Meyer

Van Daele

, et al. Exponentially fitted Runge-Kutta methods. J Comput Appl Math 2000; 125: 107–115.

Calvo

Franco

Montijano

JI.

Structure preservation of exponentially fitted Runge-Kutta methods. J Comput Appl Math 2008; 218: 421–434.

Simos

TE.

An exponentially-fitted Runge-Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions. Comput Phys Commun 1998; 115: 1–8.

Hochbruck

Lubich

Selhofer

Exponential integrators for large systems of differential equations. SIAM J Sci Comput 2006; 19: 1552–1574.

Hochbruck

Ostermann

Exponential integrators. Acta Numer 2010; 19: 209–286.

10.

De La Cruz

Biscay

Carbonell

et al. Local Linearization Runge-Kutta (LLRK) methods for solving ordinary differential equations. In: Alexandrov VN, van Albada GD, Sloot PMA and Dongarra J. (eds) Computational Science – ICCS 2006. ICCS 2006. Lecture Notes in Computer Science, vol 3991. Heidelberg, Berlin: Springer.

11.

De La Cruz

Biscay

Carbonell

, et al. A higher order local linearization method for solving ordinary differential equations. Appl Math Comput 2007; 185: 197–212.

12.

Cieśliński

Locally exact modifications of numerical schemes. Comput Math Appl 2013; 65: 1920–1938.

13.

Cieśliński

Ratkiewicz

Improving the accuracy of the discrete gradient method in the one-dimensional case. Physical Review E 2010; 81: 016704.

14.

Budd

Piggott

MD.

Geometric integration and its applications. In: Ciarlet PG (ed.), Handbook of Numerical Analysis, Vol. XI. Amsterdam: North-Holland, 2003, 11:35–139.

15.

Van Daele

Vanden Berghe

Geometric numerical integration by means of exponentially-fitted methods. Appl Numer Math 2007; 57: 415–435.

16.

Hairer

Lubich

Wanner

Geometric numerical integration illustrated by the Störmer/Verlet method. Acta Numer 2003; 12: 1–51.

17.

Chartier

Faou

Murua

An algebraic approach to invariant preserving integrators: the case of quadratic and Hamiltonian invariants. Numer Math 2006; 103: 575–590.

18.

Marsden

Lie Poisson Hamilton Jacobi theory and Lie Poisson integrators. Phys Lett A 1988; 133: 134–139.

19.

Hairer

Lubich

Wanner

Geometric numerical integration, volume 31 of springer series in computational mathematics. 2nd ed. Berlin: Springer-Verlag, 2006.

20.

Hamiltonian approach to nonlinear oscillators. Phys Lett A 2010; 374: 2312–2314.

21.

An approximate amplitude-frequency relationship for a nonlinear oscillator with discontinuity. Nonlinear Sci Lett A 2016; 7: 77–85.

22.

Ren

Gui

He's frequency formulation for nonlinear oscillators using a golden mean location. Comput Math Appl 2011; 61: 1987–1990.

23.

Sun

, et al. Active vibration control of a conical shell using piezoelectric ceramics. J Low Freq Noise V A 2017; 36: 366–375.

24.

Yildirim

Saadatnia

Askari

, et al. Higher order approximate periodic solutions for nonlinear oscillators with the Hamiltonian approach. Appl Math Lett 2011; 24: 2042–2051.

25.

Feng

Qin

Symplectic geometric algorithms for Hamiltonian systems. Heidelberg, Berlin: Springer, 2010.

26.

Leimkuhler

Reich

Simulating Hamiltonian dynamics. Cambridge: Cambridge University Press, 2004.

27.

Sanz-Serna

Calvo

MP.

Numerical Hamiltonian problems. London: Chapman & Hall/CRC, 1994.

28.

Feng

Qin

, et al. Construction of canonical difference schemes for Hamilton formalism via generating functions. J Comp Math 1989; 11: 71–96.

29.

Leok

Zhang

Discrete Hamiltonian variational integrators. IMA J Numer Anal 2011; 31: 1497–1532.

30.

Marsden

West

Discrete mechanics and variational integrators. Acta Num 2001; 10: 357–514.

31.

Beléndez

Pascual

M_Endez

, et al. Exact solution for the nonlinear pendulum. Rev Bras Ensino Fís 2007; 29: 645–648.