Periodic Solution of Nonlinear Conservative Systems

Abstract

Conservative systems represent a large number of naturally occurring and artificially designed scientific and engineering systems. A key consideration in the theory and application of nonlinear conservative systems is the solution of the governing nonlinear ordinary differential equation. This chapter surveys the recent approximate analytical schemes for the periodic solution of nonlinear conservative systems and presents a recently proposed approximate analytical algorithm called continuous piecewise linearization method (CPLM). The advantage of the CPLM over other analytical schemes is that it combines simplicity and accuracy for strong nonlinear and large-amplitude oscillations irrespective of the complexity of the nonlinear restoring force. Hence, CPLM solutions for typical nonlinear Hamiltonian systems are presented and discussed. Also, the CPLM solution for an example of a non-Hamiltonian conservative oscillator was presented. The chapter is aimed at showcasing the potential and benefits of the CPLM as a reliable and easily implementable scheme for the periodic solution of conservative systems.

1. Introduction

1.1 Hamiltonian and non-Hamiltonian conservative systems

Conservative systems can be defined as oscillating or vibrating systems in which the total energy content of the system remains constant. In order words, the total energy in the system is conserved. Ideally, such a system will continue to be in periodic oscillatory motion ad infinitum because the energy content of the system does not diminish due to the absence of dissipative force or increase due to additional energy input. However, for real cases where dissipative mechanisms such as friction or viscous damping cannot be completely eliminated, a conservative system can be thought of as one in which the energy dissipated is negligible during the time range under consideration. For example, the first few seconds of the oscillation of a simple pendulum may be considered conservative since the effect of air friction is negligible, but in the long run, the initial energy content is gradually dissipated until the pendulum comes to a halt. Other examples of practical conservative systems include mass-spring oscillator, structural elements (i.e., beams, plates, and shells), slider-crank mechanism [1], human eardrum [2], relativistic oscillator [3], planetary orbits around the sun [3], and current-carrying conductor in the electric field of an infinite rod [4]. Hence, a large number of oscillating physical systems can be studied as conservative systems.

At any point in time, the energy of a conservative system is composed of kinetic ($T$) and potential ($V$) energies except at critical points where the total energy may be only kinetic ($T_{\mathit{max}}$) or potential ($V_{\mathit{max}}$). Generally, it is expected that $T=T\left(q,\dot{q}\right)$ and $V=V\left(q,\dot{q}\right)$, where $q$ is the generalized displacement. Naturally, $q$ and $\dot{q}$ are not expected to form a product in the function $T\left(q,\dot{q}\right)$, but in some cases, they do$.$ Therefore, two types of conservative systems are distinguished namely: natural and nonnatural conservative systems. The natural conservative systems are those in which the kinetic energy can be expressed as a pure quadratic function of velocity, i.e., does not contain a product of the velocity and displacement. They are also known as Hamiltonian systems because they admit a Hamiltonian function ($H\left(q,\dot{q}\right)=T\left(q,\dot{q}\right)+V\left(q,\dot{q}\right)$) that is always constant at any point in time. While this definition of Hamiltonian systems is a physical one, a mathematical definition has been discussed by Jordan and Smith [3]. Examples of Hamiltonian systems include mass-spring oscillator, simple pendulum, and a mass attached to the mid-point of an elastic spring. On the other hand, there are conservative systems in which the kinetic energy cannot be expressed as a pure quadratic function of the velocity because the kinetic energy expression contains a product of velocity and displacement. This second group of conservative systems is referred to as nonnatural because their kinetic energy is not a pure quadratic function of velocity. Although the total energy in such systems is conserved, their Hamiltonian function ($H\left(q,\dot{q}\right)$) is not constant [4]. Hence, the nonnatural conservative systems may be referred to as non-Hamiltonian conservative systems. Examples of this category of conservative systems abound in artificial systems and include slider-crank mechanism [1], particle sliding on a vertical rotating parabola [4], pendulum attached to massless rolling wheel [4], rigid rod rocking on a circular surface without slip [4], and circular sector oscillator [5]. An important quality of the non-Hamiltonian conservative systems is that their vibration equation, which is normally derived by the Lagrangian approach, does not conform to the standard representation of conservative systems that clearly shows the restoring force. Rather, the derived vibration equation has a quadratic velocity term, which represents a coordinate-dependent parameter rather than a dissipative parameter.

1.2 Recent advances in solution schemes for nonlinear conservative oscillators

Exact analytical solutions for the nonlinear vibration models of conservative systems can be derived only in very few situations, and the solutions are usually derived in terms of special functions. Alternatively, highly accurate numerical solutions can be obtained for the nonlinear vibration model of any conservative system. However, as it is well recognized among the nonlinear science community, numerical solutions often have the limitations of lack of physical insight and convergence issues. Furthermore, there is the possibility of obtaining inaccurate convergent solutions for a nonlinear ordinary differential equation (ODE) [6], thus necessitating the independent verification of the convergent numerical solution by another numerical or analytical method. These limitations have driven the search for approximate analytical schemes capable of providing periodic solutions to nonlinear conservative oscillators. It can be rightly concluded that this search has been very fruitful considering the many approximate analytical schemes that now appear in the nonlinear science literature. The purpose of this section is to provide a brief survey of some of the notable achievements in the development of approximate analytical schemes for the periodic solution of nonlinear conservative oscillators. It should be noted that an approximate analytic method for nonlinear oscillators is considered adequate if it gives accurate predictions for the frequency-amplitude response and the oscillation history as well [7, 8].

Approximate analytical techniques to solve the nonlinear ODE governing the oscillations of a conservative system have been formulated for at least 100 years and can be classified as perturbation and nonperturbation methods. The first attempts were based on perturbation theory and are referred to as classical perturbation methods. The perturbation methods are formulated based on the concept that an unknown nonlinear system can be studied by introducing a small disturbance to a known linear system in equilibrium. For this reason, the classical perturbation methods (see Nayfeh [9] for a comprehensive treatment of classical perturbation methods) depend on the assumption of a small parameter. The problem with the small parameter assumption is that it has a small range of validity and only produces reliable solutions for cases of small-amplitude oscillations and weak nonlinearity. Nevertheless, the classical perturbation methods are still very relevant today for introducing and investigating various nonlinear concepts.

More recently, in the last four decades, a number of approximate analytical schemes have been proposed. Most of these recent schemes are nonperturbation methods, but some recent perturbation methods that attempt to improve on their classical counterparts have been formulated too. The recent perturbation methods include $\delta$-method [10], Homotopy perturbation method [11] and its variants [12, 13, 14, 15, 16, 17], modified Lindstedt-Poincare method [18, 19, 20, 21], book-keeping parameter method [22], iteration perturbation method [23], parameterized perturbation method [24], perturbation incremental method [25], and linearized perturbation method [26]. A review article on some of the recent perturbation methods has been published by He [27]. The main point of the recent perturbation methods is to deal with the issue of the small parameter in order to formulate solutions that are applicable to small- and large-amplitude oscillations and also weak and strong nonlinear oscillations. Although the higher order approximations of the recent perturbation methods have been very successful in producing accurate estimates of the frequency-amplitude response, the same cannot be said of their estimation of the oscillation history. Studies [7, 28] have shown that the higher order approximations of the recent perturbation methods produce large unbounded errors in the oscillation history during large-amplitude oscillations and are, therefore, not better than the classical perturbation methods in this regard. A plausible explanation for this observation is that it occurs because perturbation methods are based on asymptotic series that are inherently divergent for amplitudes greater than unity [28]. Therefore, it may not be possible to formulate perturbation schemes that would correctly predict the oscillation history of large-amplitude vibrations.

In contrast to the perturbation methods, the nonperturbation methods do not use any small or artificial parameter. Examples include Adomian decomposition method [29], Homotopy analysis method [30], Variational iteration method [31], Energy balance method [2] and its modifications [32, 33, 34], He Chengtian’s interpolation method [27] also called max-min approach, amplitude-frequency formulation [35], Hamiltonian approach [36], global error minimization method [37], Harmonic balance method [4] and its modifications [38, 39, 40, 41, 42], cubication methods [43, 44, 45, 46], variational methods [47, 48, 49], differential transform method [50], and continuous piecewise linearization method [8]. Nonperturbation methods also have various limitations. For instance, a study [51] showed that the Adomian decomposition method does not converge to the correct solution in some cases, and the study proposed an optimal convergence acceleration parameter to deal with this issue. Also, methods that rely on a simple harmonic approximation of the oscillation history, such as the energy balance method, amplitude-frequency formulation, Hamiltonian approach, max-min approach, and variational methods, can only give reliable estimate of the frequency-amplitude response. Sometimes these methods perform poorly in predicting the oscillation history during large-amplitude and/or strong nonlinear vibrations. Other methods that usually require high-order approximations, such as Adomian decomposition method, harmonic balance method, and variational iteration method, present algebraic complexities in their determination of higher order solutions and may be impractical for oscillators with highly complex nonlinearities such as the slider-crank mechanism [1] and the bifilar pendulum [52]. Furthermore, it has been observed that higher order estimates do not always improve the solution of the oscillation history [27]. Finally, some nonperturbation methods are heuristic in nature (e.g., energy balance method and variational methods) and require experience to choose the initial trial function and the condition for error minimization [8].

The continuous piecewise linearization method (CPLM) is an iterative analytic algorithm that was formulated to overcome most of the above challenges by providing simple and accurate solutions for the oscillation history and frequency-amplitude response of Duffing-type oscillators. In another study [53], the CPLM was modified in order to generalize it so that it can handle more complicated nonlinear conservative oscillators. Interestingly, the CPLM does not require higher order approximations or any small, artificial, or embedded parameter. Also, the algorithm is inherently stable, straightforward, and based on closed-form analytical approximations. This chapter is aimed at presenting the generalized CPLM algorithm as a veritable approach for accurate periodic solution of Hamiltonian and non-Hamiltonian conservative oscillators with complex nonlinearity. As is shown later, the CPLM retains the same simplicity in its implementation irrespective of the complexity of the nonlinear conservative system.

2. Continuous piecewise linearization method

2.1 Concept of the continuous piecewise linearization method

The main idea of the CPLM is based on the piecewise linearization of the nonlinear restoring force of a conservative oscillator. The linearization technique used by the CPLM was first applied in another algorithm [6, 54] for the solution of half-space impact models called force indentation linearization method (FILM). The FILM has been applied to formulate theoretical solutions for rigid body motions and local compliance response during nonlinear elastoplastic impact of dissimilar spheres [55]. However, because the FILM is limited to impact excitations that are nonoscillatory, it cannot be applied to nonlinear conservative oscillators. Hence, the CPLM applies the piecewise linearization technique of the FILM to provide a periodic solution for nonlinear conservative oscillators.

Essentially, the linearization technique of the CPLM involves $n$ equal discretization of the nonlinear restoring force with respect to displacement (Figure 1) and formulating a linear restoring force for each discretization. Therefore, a linear ODE can be derived for each discretization. The solution of the linear ODE approximates the solution of the original nonlinear oscillator for a time-range that is automatically determined by the CPLM and updated continuously from one discretization to the next. Details on the discretization and linearization technique of the CPLM can be found in the following references [6, 8, 54], while the applications of the CPLM to nonlinear conservative systems are presented in [8, 56].

Figure 1.

Discretization of the restoring force of a typical nonlinear oscillator.

2.2 Mathematical formulation of the continuous piecewise linearization method

The standard form for representation of a nonlinear conservative oscillator moving in the $u$-direction is given as:

$$\ddot{u}+f\left(u\right)=0,$$E1

where $f\left(u\right)$ is the nonlinear restoring force as shown in Figure 1. In Figure 1, the numbering on the horizontal axis represents the boundary points of the discretization. The $s^{\mathit{th}}$ discretization represents a general discretization with start point at $r$ and endpoint at $s=r+1$. Figure 1 shows that for each discretization, the slope of the linear approximation of the restoring force can either be positive or negative. To account for this possibility, the linearized force for the $s^{\mathit{th}}$ discretization can be expressed as:

$$F_{\mathit{rs}}\left(u\right)=\pm \left|K_{\mathit{rs}}\right|\left(u-u_r\right)+F_r,$$E2

where $K_{\mathit{rs}}=\left[f\left(u_s\right)-f\left(u_r\right)\right]/\left(u_s-u_r\right)$ is the linear slope of $F_{\mathit{rs}}\left(u\right)$ and $F_r=f\left(u_r\right)$.

Since $F_{\mathit{rs}}\left(u\right)$ is an approximation of $f\left(u\right)$ for the $s^{\mathit{th}}$ discretization, then substituting Eq. (2) in (1) gives the approximate equation of motion for each discretization as follows:

$$\ddot{u}\pm \left|K_{\mathit{rs}}\right|u=\pm \left|K_{\mathit{rs}}\right|u_r-F_r.$$E3

Eq. (10) is a nonhomogeneous linear ODE and its solution depends on whether the sign is positive or negative.

2.2.1 Solution for positive linearized stiffness

If $K_{\mathit{rs}}>0$, the solution for the displacement and velocity can be expressed as:

$$u\left(t\right)=R_{\mathit{rs}}sin\left({\omega }_{\mathit{rs}}t+{\mathrm{\Phi }}_{\mathit{rs}}\right)+C_{\mathit{rs}}$$E4a

$$\dot{u}\left(t\right)={\omega }_{\mathit{rs}}{R}_{\mathit{rs}}cos\left({\omega }_{\mathit{rs}}t+{\mathrm{\Phi }}_{\mathit{rs}}\right),$$E4b

where ${\omega }_{\mathit{rs}}=\sqrt{K_{\mathit{rs}}}$, $C_{\mathit{rs}}=u_r-F_r/K_{\mathit{rs}}$, and $R_{\mathit{rs}}={\left[{\left(u_r-C_{\mathit{rs}}\right)}^2+{\left({\dot{u}}_r/{\omega }_{\mathit{rs}}\right)}^2\right]}^{1/2}$. The initial conditions and other parameters are determined based on the oscillation stage. For the oscillation stage that moves from $+A$ to $-A$, the initial conditions for each discretization are $u_r=u_r\left(0\right)=A-r∆u$ and ${\dot{u}}_r={\dot{u}}_r\left(0\right)=-\sqrt{\left|2{\int }_A^{u_r}-f\left(u\right)\mathit{du}\right|}$, where $∆u=A/n$, and the other parameters are calculated as:

$${\mathrm{\Phi }}_{\mathit{rs}}=\left\{\begin{array}{ll}0.5\pi & {\dot{u}}_r=0\\ \pi +{tan}^{-1}\left[{\omega }_{\mathit{rs}}\left(u_r-C_{\mathit{rs}}\right)/{\dot{u}}_r\right]& {\dot{u}}_r<0\end{array}\right.$$E5a

$$∆t=\left\{\begin{array}{ll}\left(0.5\pi -{\mathrm{\Phi }}_{\mathit{rs}}\right)/{\omega }_{\mathit{rs}}& \left(u_s-C_{\mathit{rs}}\right)\ge R_{\mathit{rs}}\\ \left(0.5\pi +{cos}^{-1}\left[\left(u_s-C_{\mathit{rs}}\right)/R_{\mathit{rs}}\right]-{\mathrm{\Phi }}_{\mathit{rs}}\right)/{\omega }_{\mathit{rs}}& \left(u_s-C_{\mathit{rs}}\right)<R_{\mathit{rs}}\end{array}\right..$$E5b

For the oscillation stage that moves from $-A$ to $+A$, the initial conditions are $u_r=u_r\left(0\right)=-A+r∆u$ and ${\dot{u}}_r={\dot{u}}_r\left(0\right)=\sqrt{\left|2{\int }_A^{u_r}-f\left(u\right)\mathit{du}\right|}$; the other parameters are calculated as:

$${\mathrm{\Phi }}_{\mathit{rs}}=\left\{\begin{array}{ll}-0.5\pi & {\dot{u}}_r=0\\ {tan}^{-1}\left[{\omega }_{\mathit{rs}}\left(u_r-C_{\mathit{rs}}\right)/{\dot{u}}_r\right]& {\dot{u}}_r<0\end{array}\right.$$E6a

$$∆t=\left\{\begin{array}{ll}\left(0.5\pi -{\mathrm{\Phi }}_{\mathit{rs}}\right)/{\omega }_{\mathit{rs}}& \left(u_s-C_{\mathit{rs}}\right)\ge R_{\mathit{rs}}\\ \left(0.5\pi -{cos}^{-1}\left[\left(u_s-C_{\mathit{rs}}\right)/R_{\mathit{rs}}\right]-{\mathrm{\Phi }}_{\mathit{rs}}\right)/{\omega }_{\mathit{rs}}& \left(u_s-C_{\mathit{rs}}\right)<R_{\mathit{rs}}\end{array}\right..$$E6b

The time at the end of each discretization is $t_s=t_r+∆t$, and the end conditions $u_s$ and ${\dot{u}}_s$ are calculated by replacing $r$ with $s$ in the formulae for the initial conditions.

2.2.2 Solution for negative linearized stiffness

If $K_{\mathit{rs}}<0$, the solution for the displacement and velocity can be expressed as follows:

$$u\left(t\right)=A_{\mathit{rs}}{e}^{{\omega }_{\mathit{rs}}t}+B_{\mathit{rs}}{e}^{-{\omega }_{\mathit{rs}}t}+C_{\mathit{rs}}$$E7a

$$\dot{u}\left(t\right)={\omega }_{\mathit{rs}}\left(A_{\mathit{rs}}{e}^{{\omega }_{\mathit{rs}}t}-B_{\mathit{rs}}{e}^{-{\omega }_{\mathit{rs}}t}\right),$$E7b

where ${\omega }_{\mathit{rs}}=\sqrt{\left|K_{\mathit{rs}}\right|}$; $C_{\mathit{rs}}=u_r+F_r/\left|K_{\mathit{rs}}\right|$. Applying the initial conditions to Eqs. (7a) and (7b) gives: $A_{\mathit{rs}}=\frac{1}{2}\left(u_r+{\dot{u}}_r/{\omega }_{\mathit{rs}}-C_{\mathit{rs}}\right)$; $B_{\mathit{rs}}=\frac{1}{2}\left(u_r-{\dot{u}}_r/{\omega }_{\mathit{rs}}-C_{\mathit{rs}}\right)$. The initial and end conditions are determined in the same way as for $K_{\mathit{rs}}>0$ above. The end conditions are applied in Eq. (7a) to get the time interval for each discretization as:

$$∆t=\left\{\begin{array}{ll}\frac{1}{{\omega }_{\mathit{rs}}}{log}_e\left[\frac{\left(u_s-C_{\mathit{rs}}\right)\pm \sqrt{{\left(u_s-C_{\mathit{rs}}\right)}^2-4A_{\mathit{rs}}{B}_{\mathit{rs}}}}{2A_{\mathit{rs}}}\right]& \left(u_s-C_{\mathit{rs}}\right)>2\sqrt{A_{\mathit{rs}}{B}_{\mathit{rs}}}\\ \frac{1}{{\omega }_{\mathit{rs}}}{log}_e\left(\frac{u_s-C_{\mathit{rs}}}{2A_{\mathit{rs}}}\right)& \left(u_s-C_{\mathit{rs}}\right)\le 2\sqrt{A_{\mathit{rs}}{B}_{\mathit{rs}}}\end{array}\right..$$E8

The sign before the square root in Eq. (8) is negative for the oscillation stage that moves from $+A$ to $-A$ and vice versa. We note that if ${\dot{u}}_r=0$, then $A_{\mathit{rs}}=B_{\mathit{rs}}=\frac{1}{2}\left(u_r-C_{\mathit{rs}}\right)$ and

$$u\left(t\right)=\left(u_r-C_{\mathit{rs}}\right)cosh\left({\omega }_{\mathit{rs}}t\right)+C_{\mathit{rs}}.$$E9

Therefore,

$$∆t=\frac{1}{{\omega }_{\mathit{rs}}}{cosh}^{-1}\left(\frac{u_s-C_{\mathit{rs}}}{u_r-C_{\mathit{rs}}}\right).$$E10

2.2.3 Solution for zero linearized stiffness

In very rare situations, we may have $K_{\mathit{rs}}=0$ for one or two discretization around the turning points or relatively flat regions of the restoring force. This is likely when $∆u$ is very small, i.e., for very large $n$, and can be eliminated by increasing or decreasing $n$ slightly. However, if we want to account for $K_{\mathit{rs}}=0$, then we get [53]:

$$u\left(t\right)=H_{\mathit{rs}}+G_{\mathit{rs}}t-\frac{1}{2}{F}_r{t}^2,$$E11

where $G_{\mathit{rs}}={\dot{u}}_r+F_r{t}_r$ and $H_{\mathit{rs}}=u_r-{\dot{u}}_r{t}_r-\frac{1}{2}{F}_r{t}_r^2$. Hence, the time interval is derived from Eq. (11) as:

$$∆t=\frac{G_{\mathit{rs}}+\sqrt{G_{\mathit{rs}}^2+2F_r\left(H_{\mathit{rs}}-u_s\right)}}{F_r}.$$E12

2.3 Remarks on the CPLM algorithm

1. From the above presentation of the CPLM formulation, it is obvious that the CPLM algorithm is simple and can be implemented by undergrads without difficulty.

2. The CPLM is inherently stable and does not have convergence issues [8].

3. For few discretization, say $n\le 10$, the CPLM algorithm can be implemented with reasonable accuracy using a pocket calculator. However, the CPLM is better executed using a simple code in any programming language such as MATLAB and Mathematica or using a customized MS Excel spreadsheet.

4. When dealing with conservative oscillators with odd nonlinearity, which are symmetrical about the origin, discretization of the restoring force is only required for 0 to $A$. This means that there will be $2n$ discretizations from $-A$ to $A$.

5. The CPLM algorithm retains the same simplicity in implementation irrespective of the complexity of the restoring force. Only the $K_{\mathit{rs}}$ constant and the integral of the restoring force are to be evaluated anew for any oscillator.

6. The CPLM relies on the explicit expression of restoring force, which means that the model for the oscillator must be expressed in the form of Eq. (1). For Hamiltonian systems, the oscillator model is formulated naturally in the form of Eq. (1). For non-Hamiltonian conservative systems, the oscillator model is not formulated naturally in the form of Eq. (1). Therefore, there is a need to transform the model of non-Hamiltonian conservative systems into the form of Eq. (1) before applying the CPLM algorithm. Fortunately, this transformation only requires a simple algebraic manipulation as demonstrated in Section 4.

7. The phase equation gives the relationship between the state variables (displacement and velocity) and can be derived exactly in closed-form for all conservative systems. The CPLM extends this bilateral relationship into a tripartite one by finding the value of the corresponding independent variable (i.e., time) that matches the values of the state variables in each discretization.

8. For few discretization, say $n\le 20$, it would be necessary to extract values within each discretization in order to obtain a smooth plot of the oscillation history. The values can be extracted using the approximate solution of the displacement. However, for many discretizations, say $n\ge 50$, there is no need to extract values from any discretization.

9. The usual initial conditions investigated for nonlinear conservative oscillators are nonzero displacement and zero velocity. However, the CPLM algorithm can comfortably handle nonzero initial conditions for displacement and velocity.

3. Periodic solution of typical Hamiltonian systems

3.1 Nonlinear simple pendulum

The simple pendulum is arguably the most investigated physical system and provides very interesting insights into nonlinear phenomena. Butikov [57] calls it “an antique but evergreen physical model.” The undamped oscillation of a simple pendulum is a Hamiltonian system governed by the well-known nonlinear ODE as shown:

$$\ddot{u}+{\mathrm{\Omega }}_0^{2}sinu=0,$$E13

where $u$ is the angular displacement, ${\mathrm{\Omega }}_0=\sqrt{g/l}$, $l$ is the length of the pendulum, and $g=9.8m/s^2$. The initial conditions are given as $u\left(0\right)=A$ and $\dot{u}\left(0\right)=0$. The same initial conditions are applicable to all other oscillators discussed subsequently. The exact solution to Eq. (13) is expressed in terms of elliptic functions. The displacement and natural frequency are given as [58]:

$$u_{\mathit{ex}}\left(t\right)=2{sin}^{-1}\left[k\mathrm{sn}\left({\mathrm{\Omega }}_{0}t+K\left(k^2\right);k\right)\right]$$E14

$${\omega }_{\mathit{ex}}=\frac{\pi {\mathrm{\Omega }}_0}{2K\left(k^2\right)},$$E15

Where $\mathrm{sn}$ is the Jacobi elliptic sine function, $k=sin\left(A/2\right)$, and $K\left(k^2\right)$ is the complete elliptic integral of the first kind given as:

$$K\left(k^2\right)={\int }_0^{\pi /2}\frac{1}{\sqrt{1-k^2{sin}^2\varphi }}\mathit{d\varphi }.$$E16

From Eq. (13), the restoring force for the pendulum is $f\left(u\right)={\mathrm{\Omega }}_0^{2}sinu$ and looks like the plot in Figure 1. This means that $K_{\mathit{rs}}={\mathrm{\Omega }}_0^2\left(sin{u}_s-sin{u}_r\right)/\left(u_s-u_r\right)$ and $\dot{u}=\pm {\mathrm{\Omega }}_0\sqrt{2\left(cosu-cosA\right)}$. The initial and final velocities for each discretization are determined using the expression for $\dot{u}$, and all the other constants for the solution of each discretization are determined based on $K_{\mathit{rs}}$. A plot of the frequency-amplitude response for the simple pendulum when ${\mathrm{\Omega }}_0=1.0$ is given in Figure 2a and the corresponding error of the CPLM solution in comparison with the exact solution (Eq. (15)) is shown in Figure 2b. We see that for $A\le 178°$, the maximum error in the CPLM estimate is less than 0.45% for $n=50$ and 0.14% for $n=100$. Also, the oscillation history for moderate-amplitude ($A=45°$) and large-amplitude ($A=135°$) is shown in Figure 3 and an excellent agreement between the CPLM solution for $n=100$ and the exact solution (Eq. (14)) is observed. We noted that trigonometric nonlinearity is usually difficult to deal with and that is why the CPLM shows a relatively slow convergence to accurate results. Hence, many discretizations (e.g., $n=50-100$) are required to get an accuracy that is within 1.0% of the exact solution during large-amplitude oscillations (90$°<A<180°$) of the simple pendulum.

Figure 2.

(a) Frequency-amplitude response of the simple pendulum for $\mathit{0}°<\mathit{A}<\mathit{180}°$. (b) CPLM error analysis.

Figure 3.

Oscillation history of the simple pendulum for (a) $\mathit{A}=\mathit{45}°$ and (b) $\mathit{A}=\mathit{135}°$. CPLM estimate—Lines; exact solution—Markers.

3.2 Motion of satellite equidistant from twin stars

Consider the motion of a satellite along a path that is equidistant from two identical massive stars with mutually interacting gravitational fields. If the distance between the two stars is $2d$ and the coordinate of the satellite motion is $u$, then the equation of motion of the satellite is given as [3]:

$$\ddot{u}+\frac{2\mathit{Mu}}{{\left(d^2+u^2\right)}^{3/2}}=0,$$E17

where $M$ is the mass of a star and the restoring force is $f\left(u\right)=2\mathit{Mu}/{\left(d^2+u^2\right)}^{3/2}$. Eq. (17) shows that the satellite-star interaction results in a conservative oscillation of the satellite. Figure 4 shows the nonlinear restoring force, which is an irrational force because of the bottom square root. The restoring force spikes on both sides of the vertical axis close to the origin. The spikes indicate the point when the satellite is most influenced by the mutual gravitational field of the stars. Away from the origin, the restoring force decreases gradually and approaches the horizontal axis asymptotically. This means that the satellite is far away from the stars and experiences a much smaller gravitational force. This problem was discussed qualitatively by Jordan and Smith [3] and Arnold [59], but here, the periodic solution was investigated.

Figure 4.

Restoring force for Eq. (18): $\mathit{M}={\mathit{10}}^{\mathit{5}}\left[\mathit{kg}\right]$; $\mathit{d}=\mathit{100}\left[\mathit{m}\right]$; $\mathit{A}=\mathit{500}$.

The main CPLM constant is $K_{\mathit{rs}}=2M\left[{\left(d^2+u_s^2\right)}^{-3/2}-{\left(d^2+u_r^2\right)}^{-3/2}\right]$ and the velocity was derived as: $\dot{u}=\pm 2\sqrt{M\left[{\left(d^2+u^2\right)}^{-1/2}-{\left(d^2+A^2\right)}^{-1/2}\right]}$. The periodic solutions obtained by the CPLM and exact numerical solution are shown in Figures 5 and 6. The numerical solution was obtained by solving Eq. (17) using the NDSolve function in Mathematica™. The NDSolve function is a Mathematica subroutine for solving ordinary, partial, and algebraic differential equations numerically. In its basic form, it automatically selects the numerical method to use from a list of standard methods such as implicit Runge–Kutta, explicit Runge–Kutta, symplectic partitioned Runge–Kutta, predictor–corrector Adams, and backward difference methods. In some cases, the NDSolve function can combine two or more methods to obtain the required solution. This basic form is preferable because the NDSolve function uses the method(s) that best solves the differential equation considering accuracy and solution time. Hence, the NDSolve function was used in its basic form for all numerical solutions obtained in this chapter.

The input values used for simulation are $M={10}^5\left[\mathit{kg}\right]$ and $d=1000\left[m\right]$. In contrast to the simple pendulum, the oscillation of the satellite requires less discretization for accurate results because there is no trigonometric nonlinearity. The maximum error of the CPLM estimate for the frequency-amplitude response is less than 0.55% for $n=10$ and 0.20% for $n=20$. Significantly higher accuracies can be achieved by increasing $n$, but the results show that $n=10$ gives reasonably accurate estimates.

Figure 5.

(a) Frequency-amplitude response for satellite. (b) CPLM error analysis.

Figure 6.

Oscillation history of satellite for (a) $\mathit{A}=\mathit{50}$ and (b) $\mathit{A}=\mathit{1500}$. CPLM estimate—Lines; exact solution—Markers.

On the other hand, Figure 6 shows the oscillation history of the satellite during small-amplitude ($A=50$) and large-amplitude ($A=1500$) oscillations. The former gives a simple harmonic response with a natural frequency that is independent of the amplitude and approximately equal to $\sqrt{2M/d^3}$, while the latter exhibits an anharmonic response with a natural frequency that depends strongly on the amplitude.

3.3 Mass-spring oscillator with fractional nonlinearity

An interesting oscillator that has been the subject of several studies [60, 61, 62, 63, 64, 65] is the Hamiltonian oscillator with odd fractional nonlinearity. For the purpose of the present investigation, we consider an oscillator that is characterized by a general fractional nonlinearity as follows [65]:

$$\ddot{u}+u^{1/\left(2m+1\right)}=0,$$E18

where the restoring force, $f\left(u\right)=u^{1/\left(2m+1\right)}$, has a fractional index for all $mϵ\left\{{\mathbb{Z}}^{+}\right\}$. The main CPLM constant is evaluated as $K_{\mathit{rs}}=\left[u_s^{1/\left(2m+1\right)}-u_r^{1/\left(2m+1\right)}\right]/\left(u_s-u_r\right)$ and the velocity as $\dot{u}=\pm {\left[\left(2m+1\right)/\left(m+1\right)\right]}^{1/2}\sqrt{A^{\left(2m+2\right)/\left(2m+1\right)}-u^{\left(2m+2\right)/\left(2m+1\right)}}$. The periodic solution for the case of $m=1$, i.e., $u^{1/3}$ oscillator, is shown in Figures 7 and 8. The exact frequency-amplitude response used for the verification of the CPLM solution is [63]:

$${\omega }_{\mathit{ex}}=\frac{2\mathit{\pi \Gamma }\left(5/4\right)}{\sqrt{6}\Gamma \left(3/4\right)\Gamma \left(1/2\right)A^{1/3}}=\frac{1.070451}{A^{1/3}}.$$E19

while the exact oscillation history was obtained by the numerical solution of Eq. (18) using the NDSolve function in Mathematica™. The CPLM solution demonstrates an excellent agreement with the exact solution.

Figure 7.

(a) Frequency-amplitude response for ${\mathit{u}}^{\mathit{1}/\mathit{3}}$ oscillator. (b) CPLM error analysis.

Figure 8.

Oscillation history of ${\mathit{u}}^{\mathit{1}/\mathit{3}}$ oscillator for (a) $\mathit{A}=\mathit{0.01}$ and (b) $\mathit{A}=\mathit{10}$. CPLM estimate—lines; exact solution—markers.

Figure 7 compares the CPLM estimates of the frequency-amplitude response with Eq. (19), and the maximum error of the CPLM solution is 0.22% for $n=10$ and 0.076% for $n=20$. The error in the CPLM estimate is approximately constant for all amplitudes, and the maximum error is well below 1.0% for $n=10$. The results also reveal that the frequency approaches zero as $A\to \infty$. In Figure 8, the oscillation history for small-amplitude ($A=0.01$) and large-amplitude ($A=10.0$) oscillations is shown to exhibit similar anharmonic response, which is an indication of strong nonlinearity. Therefore, it can be concluded that the $u^{1/3}$ oscillator is highly nonlinear for all amplitudes. This quality of possessing strong nonlinearity for all amplitudes is in contrast to most Hamiltonian oscillators that are linear for small amplitudes, e.g., the oscillators considered in Sections 3.1 and 3.2 above. Another Hamiltonian oscillator that possesses strong nonlinearity for all amplitudes is the geometrically nonlinear crank [1].

4. Periodic solution of non-Hamiltonian conservative systems

The non-Hamiltonian conservative systems are generally more complex and difficult to solve compared with the Hamiltonian systems. In order to demonstrate the application of the CPLM algorithm to non-Hamiltonian conservative systems, we consider the motion of a particle on a rotating parabola. This system consists of a frictionless mass sliding on a vertical parabolic wire described by $y=q{u}^2$ for $q>0$ and rotating at a constant speed, $\mathrm{\Omega }$, about the $y$-axis (Figure 9a). The $u$-axis represents the perpendicular displacement of the mass from the $y$-axis. The kinetic and potential energies of the system are given as [4]:

$$T=\frac{1}{2}m\left[\left(1+4q^2{u}^2\right){\dot{u}}^2+{\mathrm{\Omega }}^2{u}^2\right];V=\mathit{mgq}{u}^2.$$E20

Hence, the Lagrangian is:

$$L=T-V=\frac{1}{2}m\left[\left(1+4q^2{u}^2\right){\dot{u}}^2+{\mathrm{\Omega }}^2{u}^2\right]-\mathit{mgq}{u}^2.$$E21

Figure 9.

(a) Schematic of particle on a rotating parabola. (b) Restoring force when $\mathit{q}=1.0$; $\mathit{\Lambda }=\mathit{10.0}$; $\mathit{A}=\mathit{2.0}$.

Next, we substitute Eq. (21) into the Euler–Lagrange equation to derive the equation of motion. The Euler–Lagrange equation can be written as:

$$\frac{d}{\mathit{dt}}\left(\frac{\partial L}{\partial \dot{u}}\right)-\frac{\partial L}{\partial u}=0.$$E22

Therefore, using Eqs. (21) and (22), the motion of a particle on a rotating parabola is governed by:

$$\left(1+4q^2{u}^2\right)\ddot{u}+4q^{2}u{\dot{u}}^2+\mathit{\Lambda u}=0,$$E23

where $\Lambda =2\mathit{gq}-{\mathrm{\Omega }}^2$ and the initial conditions are: $u\left(0\right)=A$ and $\dot{u}\left(0\right)=0$.

To solve Eq. (23) using the CPLM, it must be recast in the form of Eq. (1). The conservation of energy for Eq. (23) is given as [4]:

$$\left(1+4q^2{u}^2\right){\dot{u}}^2+\Lambda u^2=h,$$E24

where $h$ is a constant representing the total energy in the system. Eq. (24) confirms that the Hamiltonian, $H=T+V=\frac{1}{2}m\left[\left(1+4q^2{u}^2\right){\dot{u}}^2+\left(\Lambda +2{\mathrm{\Omega }}^2\right)u^2\right]=\frac{1}{2}m\left(h+2{\mathrm{\Omega }}^2{u}^2\right),$ is not constant. Applying the initial conditions, we get $h=\Lambda A^2$ so that ${\dot{u}}^2=\Lambda \left(A^2-u^2\right)/\left(1+4q^2{u}^2\right)$. Substituting this expression for ${\dot{u}}^2$ in Eq. (24) and after algebraic simplification, we get:

$$\ddot{u}+\frac{\Lambda \left(1+4q^2{A}^2\right)u}{{\left(1+4q^2{u}^2\right)}^2}=0.$$E25

Therefore, the restoring force is $f\left(u\right)=\Lambda \left(1+4q^2{A}^2\right)u/{\left(1+4q^2{u}^2\right)}^2$. Figure 9b shows that $f\left(u\right)$ is linear at small amplitudes and strongly nonlinear at large amplitudes. The main CPLM constant was calculated as $K_{\mathit{rs}}=\Lambda \left(1+4q^2{A}^2\right)\left[1/{\left(1+4q^2{u}_s^2\right)}^2-1/{\left(1+4q^2{u}_r^2\right)}^2\right]$, and the velocity was evaluated as $\dot{u}=\pm \sqrt{\Lambda \left(A^2-u^2\right)/\left(1+4q^2{u}^2\right)}$ . The exact time period for this oscillator can be derived in terms of elliptic function as follows [4]:

$$T_{\mathit{ex}}=4{\left[\left(1+4q^2{A}^2\right)/\Lambda \right]}^{1/2}E\left(k^2\right),$$E26

where $E\left(k^2\right)={\int }_0^{\pi /2}{\left(1-k^2{sin}^2\varphi \right)}^{1/2}\mathit{d\varphi }$ is the complete elliptic integral of the second kind and $k^2=4q^2{A}^2/\left(1+4q^2{A}^2\right)$. Then, the exact frequency was computed as ${\omega }_{\mathit{ex}}=2\pi /T_{\mathit{ex}}$, while the exact oscillation history was obtained by solving Eq. (25) numerically using the NDSolve function in Mathematica™.

A comparison of CPLM frequency estimate and the exact frequency is shown in Figure 10, while the oscillation histories for $A=0.50$ and $A=2.0$ are shown in Figure 11. As demonstrated in [4], periodic solutions for this system exist only for $\Lambda >0$. Hence, the simulations in Figures 10 and 11 were conducted for $\Lambda =10$ and $q=1.0$. An excellent agreement is observed between the CPLM estimates and the exact results. For $0<A\le 20$, the maximum error in the CPLM estimate of the frequency-amplitude response is 0.642% for $n=10$ and 0.101% for $n=20$, both of which are well below 1.0%. Also, the CPLM solution gives an accurate prediction of the strong anharmonic response in the oscillation history as shown in Figure 11.

Figure 10.

(a) Frequency-amplitude response for particle on a rotating parabola. (b) CPLM error analysis.

Figure 11.

Oscillation history of particle on a rotating parabola for (a) $\mathit{A}=\mathit{0.50}$ and (b) $\mathit{A}=\mathit{2.0}$. CPLM estimate—Lines; exact solution—Markers.

5. Concluding remarks

Conservative oscillators generally exhibit nonlinear response, and they form a large class of natural and artificially vibrating systems. Hence, the study of the dynamic response of nonlinear conservative systems is important for understanding many physical phenomena and the design of systems. The main challenge in the theoretical analysis of nonlinear conservative systems is that exact solutions are normally not available except for a few special cases where exact solutions are derived in terms of special functions.

To date, many approximate analytical methods have been formulated for the periodic solution of nonlinear conservative oscillators. This chapter provides a brief survey of the recent advances in the formulation of approximate analytical schemes and then introduced a recent approximate analytical algorithm called the continuous piecewise linearization method. The CPLM has been shown to overcome the challenges of solution accuracy and simplicity usually encountered in using most of the existing approximate analytical methods. The CPLM combines major desirable features of solution schemes such as inherent stability, accuracy, and simplicity. It is simple enough to be introduced at the undergraduate level and is capable of handling conservative oscillators with very complex nonlinearity. Conservative systems of broad interest were used to demonstrate the wide applicability of the CPLM algorithm. As demonstrated above, an accuracy of less than 1.0% relative error can be achieved for most oscillators using few discretizations, say $n\le 20$, except for oscillators with trigonometric nonlinearity where such accuracy is achieved with many discretizations. This chapter has been designed to stimulate interest in the use of CPLM for analyzing various types of conservative systems, especially those with complex nonlinearity.

Acknowledgments

The authors are grateful to the Vice-Chancellor of the University of Port Harcourt, Prof. Ndowa E.S. Lale and the Dean of Engineering, Prof. Ogbonna Joel, for approving the publication grant for this chapter. The OAPF was funded by the African Centre of Excellence in Oilfield Chemicals Research (ACE-FOR), University of Port Harcourt.