Nonlinear Violence in Nonlinear Oscillations at High Energy

Yair Zarmi

Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Midreshet Ben-Gurion, Israel.

DOI: 10.4236/am.2024.151007                

Abstract

This paper focuses on the characteristics of solutions of nonlinear oscillatory systems in the limit of very high oscillation energy, E; specifically, systems, in which the nonlinear driving force grows with energy much faster for x(t) close to the turning point, a(E), than at any position, x(t), that is not too close to a(E). This behavior dominates important aspects of the solutions. It will be called “nonlinear violence”. In the vicinity of a turning point, the solution of a nonlinear oscillatory systems that is affected by nonlinear violence exhibits the characteristics of boundary-layer behavior (independently of whether the equation of motion of the system can or cannot be cast in the traditional form of a boundary-layer problem.): close to a(E), x(t) varies very rapidly over a short time interval (which vanishes for E → ∞). In traditional boundary layer systems this would be called the “inner” solution. Outside this interval, x(t) soon evolves into a moderate profile (e.g. linear in time, or constant)—the “outer” solution. In (1 + 1)-dimensional nonlinear energy-conserving oscillators, if the solution is reflection-invariant, nonlinear violence determines the characteristics of the whole solution. For large families of nonlinear oscillatory systems, as E → ∞, the solutions for x(t) tend to common, indistinguishable profiles, such as periodic saw-tooth profiles or step-functions. If such profiles are observed experimentally in high-energy oscillations, it may be difficult to decipher the dynamical equations that govern the motion. The solution of motion in a central field with a non-zero angular momentum exhibits extremely fast rotation around a turning point that is affected by nonlinear violence. This provides an example for the possibility of interesting phenomena in (1 + 2)-dimensional oscillatory systems.

Keywords

High-Energy Oscillations, Nonlinear Violence, Boundary-Layer Characteristics

1. Introduction

There are hundreds of papers, books, theses and lecture notes that deal with the analysis of linear oscillatory systems, which are affected by a small nonlinear perturbation. The pioneering works of Poincaré [1] and Lindstedt [2] , led to the development of the Methods of Normal Forms [3] - [10] , Averaging [3] [11] [12] [13] and Multiple Time Scales [14] [15] [16] [17] [18] . Typically nonlinear oscillatory systems are the first to be studied in these references. A variety of new perturbative procedures that apply to problems, in which the nonlinearity is not small relative to the linear term, have emerged over the years. Some of these are reviewed in [19] [20] [21] [22] [23] . The asymptotic behavior of conservative oscillatory systems with polynomial nonlinearities as the highest power in the polynomial tends to infinity has been studied in the literature using a variety of expansion methods [24] - [35] . Tools for the analysis of integrals of the solutions of oscillatory systems with polynomial nonlinearity, the equation of which is

$\ddot{x}\left(t\right)+cx\left(t\right){\left|x\left(t\right)\right|}^k,c>0,k>0$ (1)

have been developed in terms of the inverse of the incomplete Beta function (Ateb(h) functions) [24] [34] [35] . Of particular relevance to the present paper is the expansion in a series of saw-tooth functions to describe the dynamics of both vibro-impact and other strongly nonlinear systems [26] [27] .

In the procedures exploited in these references the amplitude, x(t), and the period, T, are approximated by a power expansion in a small parameter through a finite order. The approximation for T leads to the emergence of secular terms, which spoil the quality of the approximation to the solution if computed for many cycles of oscillation. However, secular terms can be avoided in the analysis of energy-conserving oscillatory systems in (1 + 1) and some (1 + 2)-dimensional systems, as the period of oscillations, T, is then either known in closed form (e.g., central-field motion in a plane [36] , the Duffing [37] [38] [39] [40] , quintic [41] [42] [43] [44] and relativistic harmonic [45] [46] [47] oscillators), or can be computed to any desired level of accuracy. This idea has been exploited in the analysis of the case of a polynomial nonlinearity when both the maximal oscillation amplitude, a(E), and the highest power in the polynomial, (2N + 1), both tend to infinity [48] .

This study focuses on the behavior of nonlinear oscillatory systems in the limit of extremely high oscillation energy, E (mathematically, E → ∞), when the nonlinear driving force grows with E much faster for x(t) close to the turning point, a(E), than at any position, x(t), that is not too close to a(E). In many such systems, this behavior of the driving force determines important aspects of the solution in the limit E → ∞. In some systems it even determines the full profile of the solution in that limit. Hence, such dominance of the nonlinear part of the driving force will be called “nonlinear violence”. The main message is that when nonlinear violence dominates, the solution develops boundary-layer characteristics, which emerge independently of whether the equation of motion of the oscillator can or cannot be transformed into the traditional form of a boundary-layer system—a boundary-value problem described by an ODE, in which the highest derivative is multiplied by a small parameter, e.g.:

$\epsilon \ddot{x}\left(t\right)+x\left(t\right)+g\left(x\left(t\right),\stackrel{\dot{}}{x}\left(t\right)\right)=0,\left|\epsilon \right|\ll 1,x\left(t_1\right)=x_1,x\left(t_2\right)=x_2$ (2)

A boundary-layer analysis has been applied in hundreds of papers to the study of oscillation-relaxation systems. For the study of such problems, see, e.g., [25] [49] [50] . What characterizes the solutions of boundary-layer problems is that an approximation to the solution can be found by solving Equation (2) over two parts of the range of x. One first solves Equation (2) over the range of values of x, for which the first term in the equation, $\epsilon \ddot{x}\left(t\right)$ , is negligible owing to the smallness of ε, and $\ddot{x}\left(t\right)$ being bounded. The remainder of Equation (2) then generates an approximate solution that varies slowly. That part of the solution is called the “outer solution”; it covers most of the interval of time throughout the period of oscillations. One then solves for Equation (2) in the narrow range of t, called the “boundary layer” where $\epsilon \ddot{x}\left(t\right)$ is large owing to the rapid change in x(t). The width of this range diminishes to zero as ε → 0 and $\ddot{x}\left(t\right)$ increases indefinitely as ε → 0. That part of the solution is called the “inner solution”. Both “inner” and “outer” solutions are determined up to some free parameters. The approximate solution for the system over the whole range of x is found by a “matching” procedure: requiring that the values of the “inner” and “outer” solutions and of their first derivatives coincide at some matching point, t_m. This determines the free parameters, yields the value of the matching point, t_m and the narrow range in x in the vicinity f t_m, over which x(t) varies rapidly. As ε → 0, the variation of x(t) throughout the boundary layer becomes faster and the width of the layer tends to zero.

One purpose of this paper is to show that solutions of nonlinear oscillatory systems that are dominated by nonlinear violence exhibit characteristics that are similar to those of boundary-layer problems. When nonlinear violence occurs in the vicinity of a turning point, the solution varies rapidly around that point over a short time interval, the length of which tends to zero as E → ∞. In boundary-layer terminology, this would be called the “inner solution”. Beyond this interval, the solution rapidly evolves into a moderate profile, often, well approximated by linear time dependence, in other systems—by a constant. In boundary-layer terminology, this would be called the “outer solution”. Exploiting the analysis tools developed in boundary-layer theory facilitates understanding of the characteristics of the solutions, in particular, in the limit of very high oscillation energy.

Section 2 deals with (1 + 1)-dimensional conservative systems in the limit of total energy, E → ∞. Systems not transformable to boundary-layer formulation are analyzed in Section 2.1. The first to be discussed, are systems, the solutions of which are symmetric under space inversion and tend to saw-function patterns. The general case of potential energy, V(x), that is symmetric under reflection in x, when V(x) grows sufficiently rapidly near the turning points, x = ±a(E) is studied in Section 2.1.1. An example of a Gaussian potential is presented in Section 2.1.2. The case of a polynomial driving force with constant coefficients with a positive coefficient multiplying the leading power is reviewed in Section 2.1.3. The case of the extreme relativistic limit of the relativistic harmonic oscillator is discussed in Section 2.2. Unlike all previously discussed systems, in the extreme relativistic limit, the dynamical equation of the relativistic harmonic oscillator is transformable to a boundary-layer problem. In all the system discussed, the solutions are governed by nonlinear violence. The growth characteristics of V(a(E)) lead to a profile of x(t), which, to a good approximation, is comprised of rising and decreasing segments with linear time dependence throughout most of each half of the period, so that $\stackrel{\dot{}}{x}\left(t\right)$ is very close to a constant value rapidly changing sign over the vanishingly small interval around the turning points. Around the turning points, the solution is highly curved. Consequently, As E → ∞, x(t) tends to a periodic saw-tooth profile, the velocity, $\stackrel{\dot{}}{x}\left(t\right)$ ,—to a periodic step-function, the acceleration, $\ddot{x}\left(t\right)$ ,—to a periodic sequence of δ-functions and the phase-space plot, $\stackrel{\dot{}}{x}\left(t\right)$ vs. x(t),—to a rectangle.

A different effect of nonlinear violence emerges in some systems with polynomial nonlinearities, in which the coefficients of one or more of the lower powers in the polynomial grow with the maximal amplitude, a(E), so rapidly that they compete with the effect of the leading power of the polynomial when the solution, x(t), is close to a(E). For such systems, x(t) tends to a step function when a(E) grows indefinitely. The Duffing equation, in which the coefficient of the linear term grows with energy, is discussed in Section 2.3 because of its special physical significance. It was obtained by an ingenious transformation of the relativistic harmonic oscillator [46] [47] [51] [52] [53] .

The case, in which V(x) is not symmetric under reflection in x, and nonlinear violence emerges only near one turning point is discussed in Section 3.

The motion of a particle in a central field, V(r) (r is the distance of the particle from the origin) in (1 + 2) dimensions, when E → ∞ and the particle has non-zero angular momentum, opens the door to interesting phenomena in higher space dimensions. It is discussed in Section 4. The classical case of motion under a gravitational field, V(r) µ −1/r, follows a single “oscillation” with an infinite period of motion [36] . At high E, if the system has a non-zero angular momentum, l, the dominance of the latter in the effective potential leads to nonlinear violence around the single turning point, the smallest radius, r_min. If, instead, V(r) is chosen to grow rapidly at large r, the motion is periodic, with turning points at the smallest radius, r_min, and the largest radius, r_max. For E >> 1, the solution exhibits nonlinear violence at both turning points. At r_min, this behavior is caused by the dominance of the angular-momentum term in the effective potential (as E → ∞, r_min → 0). At r_max, this behavior is caused by the dominance of the potential, V(x), (as E → ∞, r_max → ∞). The boundary-layer characteristics of r(t) around r_min, leads to a rapid jump in the angle of rotation around the origin, φ, when crossing r_min, within a short time interval, which vanishes when E → ∞. This jump is followed by very slow variation of φ with time away from r_min. Summarizing comments are presented in Section 5.

In most systems discussed, when a turning point that is affected by nonlinear violence is crossed (crossing time t₀), as E → ∞, the velocity jumps by a finite value, Δv. As a result, the acceleration tends to Δv ∙ δ(t − t₀). For each system discussed, it is shown that, within the accuracy of calculational procedure, this jump is reproduced by the “inner solution” of the boundary layer-like analysis.

2. Nonlinear Violence in (1 + 1) Dimensional Oscillatory Systems

2.1. System Not Transformable to Boundary Layer Formulation

Consider an energy-conserving oscillatory system:

$\ddot{x}\left(t\right)+V^{\prime }\left(x\left(t\right)\right)=0$ (3)

where V(x) → ∞ as x → ∞ and the maximal amplitude, a(E) tends to infinity as the total energy, E → ∞. To ensure oscillatory behavior, V(x), must be positive at both turning points. Nonlinear violence emerges in the vicinity of a turning point, x = a(E), if V(x) grows very rapidly in the vicinity of that point.

2.1.1. Reflection Invariant Potential—General Case

To see how a boundary layer-like behavior emerges when V(x) is symmetric under reflection in x, scaled coordinates will be used:

$x\left(t\right)=a\eta \left(\theta \right),\theta =\frac{2\pi t}{T}$ (4)

The motivation for the choice of the angle variable, θ, is that the motion is periodic with period T. As T can be found either in closed form, or computed to numerically to any desired accuracy (see discussion in Section 1), in the analysis presented here, T is treated as a known constant. Furthermore, as the oscillations vary in the range −a ≤ +a, to obtain a standard limiting profile form that applies to many systems, the coordinate, x(t) is replaced by the bounded scaled coordinate, η, −1 ≤ η ≤ +1, with period (2π) in θ. Equation (3) becomes a boundary value problem for η:

${\eta }^{″}\left(\theta \right)+\frac{1}{a}{\left(\frac{T}{2\pi }\right)}^2{V}^{\prime }\left(a\eta \left(\theta \right)\right)=0,\eta \left(0\right)=-\eta \left(\pi \right)=\eta \left(2\pi \right)=1$ (5)

The expression for the period, T is:

$\begin{array}{l}T=4{\displaystyle \underset{0}{\overset{a}{\int }}\left(1/\sqrt{2\left(E-V\left(x\right)\right)}\right)\text{d}x}=a\left(\sqrt{8}/\sqrt{V\left(a\right)}\right)I\\ I={\displaystyle \underset{0}{\overset{1}{\int }}\left(1/\sqrt{1-\left(V\left(a\eta \right)/V\left(a\right)\right)}\right)\text{d}\eta }\end{array}$ (6)

Often the integral, I, defined in Equation (6) is O(1). For example, if the potential function obeys the requirement:

$V\left(a\eta \right)/V\left(a\right)\underset{\begin{array}{l}a\to \infty \\ \left|\eta \right|<1,\text{fixed}\end{array}}{\to }C<1$ (7)

where C is a constant, independent of a, then

$I\underset{a\to \infty }{\to }\frac{1}{\sqrt{1-C}}$ (8)

In many physical systems, C = 0, and then I → 1. When Equation (7) holds, T vanishes in proportion to (V(a))^−1/2 = E^−1/2.

Exploiting Equation (6), Equation (5) becomes:

${\eta }^{″}\left(\theta \right)+2\frac{I^2}{{\pi }^2}\frac{a{V}^{\prime }\left(a\eta \left(\theta \right)\right)}{V\left(a\right)}=0$ (9)

When nonlinear violence dominates, Equation (9) leads to the emergence of behavior of the solution, which looks like the “inner” and “outer” solutions of boundary-layer problems. This happens if

$I^2\frac{a{V}^{\prime }\left(a\eta \right)}{V\left(a\right)}\underset{\begin{array}{l}a\to \infty \\ \left|\eta \right|<1,\text{fixed}\end{array}}{\to }0$ (10)

$I^2\frac{a{V}^{\prime }\left(a\right)}{V\left(a\right)}\underset{a\to \infty }{\to }\infty$ , (11)

Note the dramatic difference between Equations (10) and (11). Equation (10) ensures that, as long as |η| < 1 (specifically, (1 − |η|) >> O(1/a)), namely, not too close to the turning points, the second term in Equation (9) becomes negligible, so that ${\eta }^{″}\left(\theta \right)$ is very small. Thus, Equation (11) ensures that, except near the turning points, Equation (9) tends to

${\eta }^{″}\left(\theta \right)=0,\left(a\to \infty ,\left|1-\eta \right|\gg O\left(1/a\right)\right)$ . (12)

Equation (11) yields a linear profile for the “outer solution”:

${\eta }_{\text{outer}}\left(\theta \right)=A\theta +B$ . (13)

Consider the interval, 0 ≤ θ ≤ π, with initial conditions: η(0) = 1 and η'(0) = 0. Owing to symmetry and periodicity of the solution, one has η_outer(π/2) = 0. Also, although η_outer may miss the actual solution very close θ = 0 and π, up to errors of at most O(1/a), it must reach close to +1 as θ → 0 and -1 as θ → π: Hence, Equation (13) must be:

${\eta }_{\text{outer}}\left(\theta \right)\approx A\left(\pi /2+\delta \left(a\right)-\theta \right),A=\left(2/\pi \right)+\epsilon \left(a\right),\delta \left(a\right),\epsilon \left(a\right)\underset{a\to \infty }{\to }0$ (14)

On the other hand, Equation (11) ensures that, very close to the turning point, ${\eta }^{″}\left(\theta \right)$ is very large, so that the maximum η(0) = +1 is very sharp. η(θ) varies rapidly over a very short range in θ. Hence, θ = 0 is the tip of a boundary layer. As η(0) = +1 is a point of a symmetric maximum for η(θ), approximate η(θ) for small θ by

${\eta }_{\text{inner}}\left(\theta \right)=1-b{\theta }^2+c{\theta }^4+O\left({\theta }^6\right)$ . (15)

Choosing the nomenclature “inner” for this part of the solution is based on the observation presented below, that it varies rapidly over a narrow range in θ in the vicinity of θ = 0. This is just the characteristic of the inner solution in a boundary-layer problem. Replace η(θ) in Equation (9) by η_inner(θ) and expand the result around θ = 0 through θ⁴ to find:

$b=\frac{I^2}{{\pi }^2}\frac{a{V}^{\prime }\left(a\right)}{V\left(a\right)},c=O\left({\left(I^2\frac{a{V}^{\prime }\left(a\right)}{V\left(a\right)}\right)}^2\right)$ . (16)

Requiring that the inner and outer solutions coincide and so do their derivatives at a matching point, θ_m, one obtains for θ_m, and for A (see Equation (13)):

$\begin{array}{l}{\theta }_{\text{m}}=\frac{\pi }{I^2}\frac{V\left(a\right)}{a{V}^{\prime }\left(a\right)}+O\left({\left(\frac{V\left(a\right)}{a{I}^2{V}^{\prime }\left(a\right)}\right)}^2\right)\\ A=\frac{2}{\pi }\left(1+\frac{1}{I^2}\frac{V\left(a\right)}{a{V}^{\prime }\left(a\right)}+O\left({\left(\frac{1}{I^2}\frac{V\left(a\right)}{a{V}^{\prime }\left(a\right)}\right)}^2\right)\right)\end{array}$ . (17)

Equation (11) ensures that the range in θ, over which the “inner solution” evolves, is small and vanishes for a → ∞. Consequently, the matching point, θ_m of Equation (17), shrinks to zero and, concurrently, that A tends to (2/π) as a → ∞. With Equations (10) and (11) satisfied, as a → ∞, η(θ) tends to a periodic saw-tooth profile: Linear segments connect the turning points (η(2nπ) = +1, η((2n + 1)π) = −1), the velocity, η'(θ), tends to a periodic step-function profile, alternating between ±(2/π) and the phase-space plot, η'(θ) vs. η(θ), tends to a rectangle.

δ-function profile for acceleration The periodic step-function limit of the profile of η'(θ) as a → ∞ implies that, η''(θ), tends to a periodic sequence of δ-functions:

${\eta }^{″}\left(\theta \right)\underset{a\to \infty }{\to }{\left(-1\right)}^{n+1}\frac{4}{\pi }\delta \left(\theta -n\pi \right)$ . (18)

The δ-functions arise from the cumulative change in velocity, η', across a turning point. Within a very a good approximation, the contribution of the “inner solution” generates the jump in η' as the turning point is crossed, as a → ∞. Using η_inner(θ) of Equation (14) through O(θ²) for η(θ) as well as Equations (15) and (16), yields the total change in the slope:

${\eta }^{\prime }\left({\theta }_{\text{m}}\right)-{\eta }^{\prime }\left(-{\theta }_{\text{m}}\right)={\displaystyle \underset{-{\theta }_{\text{m}}}{\overset{+{\theta }_{\text{m}}}{\int }}{\eta }^{″}\left(\theta \right)\text{d}\theta }\approx {\displaystyle \underset{-{\theta }_{\text{m}}}{\overset{+{\theta }_{\text{m}}}{\int }}{{\eta }^{″}}_{\text{inner}}\left(\theta \right)\text{d}\theta }=-2b{\theta }_{\text{m}}=-\frac{4}{\pi }$ . (19)

Note that the result of Equation (19) does not depend on the asymptotic behavior of the period, T or the integral I (see Equation (6) and the following discussion).

2.1.2. Gaussian Potential

As an example to the statements made in Section 2.1.1, consider a Gaussian potential,

$V\left(x\right)=\frac{1}{2}{\text{e}}^{x^2}$ (20)

Equations (10) and (11) are satisfied for this potential as a → ∞.

Using the scaled coordinates of Equation (4), the period is given by:

$T=4a{\text{e}}^{-a^2/2}I,I={\displaystyle \underset{0}{\overset{1}{\int }}\left(1/\sqrt{1-{\text{e}}^{-a^2\left(1-{\eta }^2\right)}}\right)\text{d}\eta }=1+{\text{e}}^{-a^2}\frac{Erf\left[ia\right]}{4ia}+\cdots \underset{a\to \infty }{\to }1$ (21)

For the error function, see [54] . As a grows, I tends to 1 rapidly. Already for a = 10, one finds I = 1.007. In addition, Equation (17) yields that the matching point, θ_m , is O(1/a²).

The equation for η(θ) is:

${\eta }^{″}\left(\theta \right)+\frac{4a^2{I}^2}{{\pi }^2}{\text{e}}^{-a^2\left(1-\eta {\left(\theta \right)}^2\right)}\eta \left(\theta \right)=0$ . (22)

Figures 1-3 present numerical solutions of Equation (22), showing, respectively, how η(θ) tends to a periodic saw-tooth profile, η'(θ)—to a periodic step-function and the phase plot, η'(θ) vs. η(θ), to a rectangle. Let us see how this happens.

Figure 1. Approach of η(θ), solution of Equation (22), to periodic saw-tooth profile as maximal amplitude, a, is increased.

Figure 2. Approach of η'(θ), solution of Equation (22), to periodic step-function profile as maximal amplitude, a, is increased.

Figure 3. Approach of η(θ) − η'(θ) phase-space plot of solution of Equation (22), to rectangle as maximal amplitude, a, is increased.

At the turning point |η| = 1 and |η'| = 0, energy conservation for Equation (22) yields:

${\eta }^{\prime }\left(\theta \right)=\pm \frac{2}{\pi }\sqrt{1-{\text{e}}^{-a^2\left(1-\eta {\left(\theta \right)}^2\right)}}I$ . (23)

Equation (23) implies that max|η'(θ)| (attained at θ = π/2 and 3π/2, where η(θ) = 0), is:

$\max \left|{\eta }^{\prime }\left(\theta \right)\right|=\frac{2}{\pi }\sqrt{1-{\text{e}}^{-a^2}}I\underset{a\to \infty }{\to }\frac{2}{\pi }$ . (24)

The limit of Equation (24) also holds for η'(θ) of Equation (23) as long as (1 − η(θ)²) = O(1):

${\eta }^{\prime }\left(\theta \right)\underset{\begin{array}{l}a\to \infty \\ 1-\eta {\left(\theta \right)}^2=O\left(1\right)\end{array}}{\to }=\pm \frac{2}{\pi }$ . (25)

Thus, except near the turning points, where η(θ) ≈ ±1, the profile of the “outer solution” tends to a straight line, which, owing to the symmetry and periodicity of the solution, must pass through zero at θ = π/2:

${\eta }_{\text{outer}}\left(\theta \right)\underset{\begin{array}{l}a\to \infty \\ 1-\eta {\left(\theta \right)}^2=O\left(1\right)\end{array}}{\to }\pm 2/\pi \left(\pi /2-\theta \right)$ . (26)

Only near the turning points, for (1 − η(θ)²) = O(1/a²), is this limit not attained, and the “inner solution” evolves. For η(0) = 1 and η'(0) = 0 at the turning point at θ = 0, the solution is maximal and symmetric in deviations around θ = 0. Hence, expand η(θ):

${\eta }_{\text{inner}}\left(\theta \right)=1-b{\theta }^2+O\left({\theta }^4\right)$ . (27)

Substituting Equation (27) in Equation (22), and expanding the equation through the O(θ²) yields:

$b=\frac{{\text{e}}^{a^2}{T}^2}{8{\pi }^2}=\frac{2}{{\pi }^2}{a}^2{I}^2$ . (28)

${{\eta }^{\prime }}_{\text{inner}}\left(\theta \right)=-\frac{4a^2{I}^2}{{\pi }^2}\theta +O\left({\theta }^3\right)$ . (29)

Using Equations (21) and (27)-(29) yields the “inner solution”, whereas the “outer solution” is determined up to an overall coefficient:

${\eta }_{\text{inner}}=1-\frac{2a^2{I}^2}{{\pi }^2}{\theta }^2+O\left({\theta }^4\right),{\eta }_{\text{outer}}=A\left(\frac{\pi }{2}-\theta \right)$ . (30)

Requiring that the two solutions and their derivatives coincide at a matching angle, θ_m, in the vicinity of θ = 0, one obtains:

${\theta }_{\text{m}}\underset{a\gg 1}{\approx }\frac{\pi }{2a^2{I}^2}\left(1+O\left(1/a^2\right)\right)$ , (31)

$A=\frac{2}{\pi }\left(1+O\left(1/a^2\right)\right)$ . (32)

θ_m, the range, over which the “inner solution” evolves, shrinks to zero and A tends to (2/π) at an O(1/a²) rate.

In summary, as a → ∞, the solution for η(θ) tends to a periodic saw-tooth profile, which varies linearly from +1 at θ = 0, to −1 at θ = π and back to +1 at θ = 2π. η'(θ) tends to a periodic step function, with values of −(2/π) over 0 ≤ θ ≤ π and +(2/π) over π ≤ θ ≤ 2π.

δ-function profile for acceleration As in Section 2.1.1, η''(θ) tends to a periodic sequence of δ-functions with alternating signs (Equation (18)). Again, the jump in η'(θ) as the turning point is crossed is well approximated by the cumulative change of η'(θ) over the “inner solution”. Substituting η_inner(θ) of (24)-(29) for η(θ) through O(θ²) one obtains:

${\eta }^{\prime }\left({\theta }_{\text{m}}\right)-{\eta }^{\prime }\left(-{\theta }_{\text{m}}\right)\approx {\displaystyle \underset{-{\theta }_{\text{m}}}{\overset{+{\theta }_{\text{m}}}{\int }}{{\eta }^{″}}_{\text{inner}}\left(\theta \right)\text{d}\theta }=-2b{\theta }_{\text{m}}=-\frac{4}{\pi }$ . (33)

This is the total change of the slope of the “outer solution” as the turning point is crossed.

2.1.3. Polynomial Nonlinearity: Constant Coefficients

This Section focuses on the case of a polynomial V(x):

$\ddot{x}\left(t\right)+x\left(t\right)+{\displaystyle \underset{n=1}{\overset{N}{\sum }}{a}_{2n+1}x{\left(t\right)}^{2n+1}}=0,a_{2N+1}>0$ . (34)

When all coefficients, a_2n+1, in Equation (34), are constant, the leading power, x^2N+1, dominates around the turning points in the limit of very high energy, E, as it grows like a^2N+1. Still, the nonlinearity in Equation (34) does not obey conditions (10), (11) and (16) for fixed (2N + 1). The asymptotic saw-tooth profile of the solution emerges only in the double limit E → ∞ and N → ∞. As this case has been discussed in detail in [34] , the results are briefly reviewed here. Defining a small parameter,

$\mu =1/\left({\left(a_{2N+1}\right)}^{1/2N}a\right)\ll 1$ , (35)

the period of oscillation, T, is found to be O(μ^N):

$T={\mu }^{N}G\left(N,\mu \right),G=4\sqrt{N+1}\left(1+O\left(1/N\right)\right)+O\left({\mu }^{2p}\right)$ , (36)

where (2p + 1) is the power of the next-to the leading term.

Under the scaling of Equation (4)), Equation (34) becomes:

${\eta }^{″}\left(\theta \right)+{\left(\frac{G}{2\pi }\right)}^2\left\{{\mu }^{2N}\eta \left(\theta \right)+{\displaystyle \underset{n=1}{\overset{N}{\sum }}\frac{a_{2n+1}}{{\left(a_{2N+1}\right)}^{n/N}}{\mu }^{2\left(N-n\right)\mu }\eta {\left(\theta \right)}^{2n+1}}\right\}=0$ . (37)

In the double limit of a → ∞ (μ → 0) and N → ∞, for all |η(θ)| < 1, Equation (37) is reduced to

$a\to \infty :{\eta }^{″}\left(\theta \right)+{\left(G/\left(2\pi \right)\right)}^2\eta {\left(\theta \right)}^{2N+1}=0\underset{\begin{array}{l}N\to \infty \\ \left|\eta \left(\theta \right)<1\right|\end{array}}{\Rightarrow }{\eta }^{″}\left(\theta \right)=0$ , (38)

leading to the emergence of the “outer solution”.

The boundary-layer-like analysis yields:

$\begin{array}{l}{\eta }_{\text{outer}}\left(\theta \right)=\left(2/\pi \right)/\left(\left(\pi /2\right)-\theta \right)\left(1+O\left(1/\left(N+1\right)\right)\right)\\ {\eta }_{\text{inner}}\left(\theta \right)=1-\frac{2\left(N+1\right)}{{\pi }^2}{\theta }^2\left(1+O\left(1/\left(N+1\right)\right)\right)\\ {\theta }_{\text{m}}\underset{a\to \infty }{\to }\frac{\pi }{2N+2}\left(1+O\left(1/\left(N+1\right)\right)\right)\end{array}$ . (39)

Nonlinear violence emerges in the double limit of a → ∞ and N → ∞, leading to the same periodic saw-tooth profile for η(θ) as in previously discussed systems.

δ-function profile for acceleration

As in Sections 2.1.1 and 2.1.2, η''(θ) tends to a periodic sequence of δ-functions (Equation (18)). Again, the jump in η'(θ) generated by the δ-function is well approximated by cumulative change of η'(θ) over “η_inner” of Equation (39):

${\eta }^{\prime }\left({\theta }_m\right)-{\eta }^{\prime }\left(-{\theta }_m\right)\approx {\displaystyle \underset{-{\theta }_m}{\overset{+{\theta }_m}{\int }}{{\eta }^{″}}_{\text{inner}}\left(\theta \right)\text{d}\theta }\underset{a\to \infty }{\to }-\frac{4}{\pi }\left(1+O\left(1/\left(N+1\right)\right)\right)\underset{N\to \infty }{\to }-\frac{4}{\pi }$ (40)

2.2. Relativistic Harmonic Oscillator: Transformable into Boundary-Layer Problem

2.2.1. Introductory Comments

In the extreme relativistic limit of the relativistic harmonic oscillator,

$m{\omega }_0^{2}x+m\frac{\ddot{x}}{{\left(1-{\left(\stackrel{\dot{}}{x}/c\right)}^2\right)}^{3/2}}=0$ . (41)

x(t) tends to a periodic saw-tooth profile. Over the years, this has been demonstrated through numerical solutions of the Equation (41) [45] [46] [47] . A closed-form expression for the solution of Equation (41) [46] [51] [52] [53] , confirms the numerical results. For c, the speed of light in vacuum, there is little interest in the extreme relativistic limit of Equation (41): Near the turning points, the acceleration is so high, that General Relativistic effects become important. Hence, to perform an experiment, the oscillator must be placed in outer space far away from any large mass. Interest in the relativistic harmonic oscillator has recently emerged in low-temperature physics, in the oscillations of ions in a Bose-Einstein condensate [55] , where the speed of light is smaller than in vacuum by orders of magnitude.

To present numerical solutions of Equation (41) it pays to use dimensionless variables:

$\tilde{t}={\omega }_{0}t,\tilde{x}\left(\tilde{t}\right)=\left({\omega }_0/c\right)x\left(t\right)$ . (42)

Equation (41) is transformed into

$\tilde{x}\left(\tilde{t}\right)+\frac{\ddot{\stackrel{˜}{x}}\left(\tilde{t}\right)}{{\left(1-{\left(\stackrel{\dot{}}{\tilde{x}}\left(\tilde{t}\right)\right)}^2\right)}^{3/2}}$ (43)

Define

$\beta =\max \left(\left|\stackrel{\dot{}}{x}\left(t\right)\right|/c\right)=\max \left(\stackrel{\dot{}}{\tilde{x}}\left(\tilde{t}\right)\right)$ . (44)

An example of phase-space plots for $\stackrel{\dot{}}{\tilde{x}}\left(\tilde{t}\right)$ vs. $\tilde{x}\left(\tilde{t}\right)$ is shown in Figure 4 for β = 0.3, 0.7, 0.9999. As β → 1 (E → ∞), the plots tend to a rectangle. Such plots have been presented in many papers (see, e.g., [31] [32] [33] [37] [38] [39] .)

That the solution of Equation (41) tends in the extreme relativistic limit to the same limiting saw-tooth profile as in previous Sections may not be surprising, as $\stackrel{\dot{}}{x}\left(t\right)$ tends to ±c, the speed of light, throughout most of the period of oscillations. The characteristic that is common to the present system and the previously discussed systems, is that the solution of Equation (41) also exhibits a boundary-layer-like structure: most of the time, the solution is well approximated by the linear time dependence of an “outer solution”. Around each turning point, the solution of x(t) has a strongly curved profile over a very small fraction of the total period, T, a fraction that tends to zero as E → ∞. The curved part of the solution is very well approximated by an “inner solution”. Dominance of the nonlinearity near the turning points ensures that the asymptotic form of the solution of Equation (41) is identical to that found in Sections 2.1.1 - 2.1.3. This happens despite the fact that there are qualitative differences between the relativistic harmonic oscillator and the previously discussed systems. In the previous systems, one has:

$E\to \infty :a\left(E\right)\to \infty ,\max \left|\stackrel{\dot{}}{x}\left(t\right)\right|\propto \sqrt{E}\to \infty ,T\to 0$ , (45)

where E is the total energy. In the case of the relativistic harmonic oscillator, one has:

$E\to \infty :a\left(E\right)\to \infty ,\max \left|x^{\prime }\left(t\right)\right|\to c,T\to \infty$ . (46)

Moreover, in the systems discussed in previous subsections, nonlinear violence emerged without resorting to the transformation of the dynamical equations into traditional boundary-layer systems of the form of Equation (2). In the case of Equation (41), nonlinear violence emerges in the extreme relativistic limit, through such a transformation.

2.2.2. Boundary Layer-Like Analysis

Energy conservation for Equation (41) yields:

$a\left(E\right)=\max \left|x\right|=\sqrt{2}\frac{c}{{\omega }_0}\frac{\sqrt{1-\sqrt{1-{\beta }^2}}}{{\left(1-{\beta }^2\right)}^{1/4}},\max \left|\stackrel{\dot{}}{x}\right|=\beta c$ . (47)

The transformation of Equation (4) yields a periodic boundary value problem:

Figure 4. Phase-space plots of solution of Equation (43). From inward to outward: β = 0.3, 0.7, 0.9999.

$\eta \left(\theta \right)+{\left(\frac{2\pi }{{\omega }_{0}T}\right)}^2\frac{{\eta }^{″}\left(\theta \right)}{{\left(1-{\left(\frac{2\pi x_{\max }}{cT}\right)}^2{\eta }^{\prime }{\left(\theta \right)}^2\right)}^{3/2}}=0,\begin{array}{l}\eta \left(0\right)=-\eta \left(\pi \right)=\eta \left(2\pi \right)=1\hfill \\ \eta \left(\pi /2\right)=0\hfill \end{array}$ . (48)

η(θ) is (2π)-periodic in θ and bounded: |η(θ)| ≤ 1. The period, T, is given by [31] :

$\begin{array}{l}T=\frac{1}{{\omega }_0}4\sqrt{2}\left(\sqrt{1+\frac{1}{\sqrt{1-{\beta }^2}}}E\left(k\right)-\left(1/\sqrt{1+\frac{1}{\sqrt{1-{\beta }^2}}}\right)K\left(k\right)\right)\\ k^2=\frac{1-\sqrt{1-{\beta }^2}}{1+\sqrt{1-{\beta }^2}}\end{array}$ , (49)

where E(k) and K(k) are the complete elliptic integrals [54] :

$E\left(k\right)={\displaystyle \underset{0}{\overset{\pi /2}{\int }}\sqrt{1-k^2\sin {\phi }^2}\text{d}\phi },K\left(k\right)={\displaystyle \underset{0}{\overset{\pi /2}{\int }}\left(1/\sqrt{1-k^2\sin {\phi }^2}\right)\text{d}\phi }$ . (50)

Applying energy conservation to Equation (48) and exploiting properties of the elliptic integrals [54] the maximal values of η and η' are found to be:

${\eta }_{\max }\equiv \max \left(\left|\eta \right|\right)=1,{{\eta }^{\prime }}_{\max }\equiv \max \left|\left({\eta }^{\prime }\left(\theta \right)\right)\right|=\frac{\beta cT}{2\pi x_{\max }}\underset{\beta \to 1}{\to }\frac{2}{\pi }$ . (51)

To measure the proximity to the extreme relativistic limit, define a small parameter, μ:

${\beta }^2=1-{\mu }^4,\left(\mu \ll 1\right)$ . (52)

Expansion in power of μ, (for elliptic integrals, see [54] ) yields to lowest order in μ:

$\begin{array}{l}\eta \left(\theta \right)+\frac{{\pi }^2}{8}{\mu }^2\left(1+O\left({\mu }^4\log \mu \right)\right)\frac{{\eta }^{″}\left(\theta \right)}{{\left(1-\frac{{\pi }^2}{4}\left(1-{\mu }^2+O\left({\mu }^4\right)\right){\eta }^{\prime }{\left(\theta \right)}^2\right)}^{3/2}}=0\\ \eta \left(0\right)=1,{\eta }^{\prime }\left(0\right)=0\end{array}$ . (53)

The system exhibits the characteristics of a boundary layer problem (see Equation (2)): the highest derivative in the Equation (53) is multiplied by a small parameter.

To find the inner solution corresponding to Equation (53) near the turning point, θ = 0, (η(0) = 1, η'(0) = 0), observe that η' is very close to zero there. Replacing η' by zero in Equation (53), the O(μ²) approximation to Equation (53) becomes:

$\eta \left(\theta \right)+\frac{{\pi }^2}{8}{\mu }^2{\eta }^{″}\left(\theta \right)=0$ . (54)

For the boundary conditions at the turning point, Equation (54) is solved by

${\eta }_{\text{inner}}\left(\theta \right)=\cos \left(\frac{2\sqrt{2}}{\pi }\frac{\theta }{\mu }\right)$ . (55)

To obtain the “outer solution” we focus on the vicinity of θ = (π/2). As the chosen turning point is at θ = 0, two observations hold near θ = (π/2): η(θ) ≈ 0 and η' ≈ ( $-{{\eta }^{\prime }}_{\max }$ ). Substituting η(θ) = 0 in Equation (53), yields the approximate equation;

${\eta }^{″}=0\Rightarrow {\eta }_{\text{outer}}={{\eta }^{\prime }}_{\max }\left(\frac{\pi }{2}-\theta \right)$ . (56)

Exploiting the leading-order term in ${{\eta }^{\prime }}_{\max }$ by expanding Equation (51) in powers of μ, yields:

${{\eta }^{\prime }}_{\max }\equiv \max \left|\left({\eta }^{\prime }\left(\theta \right)\right)\right|=\frac{2}{\pi }+\frac{{\mu }^2}{\pi }+O\left({\mu }^4\log \mu \right)$ . (57)

Using Equations (55) - (57), an O(μ²) approximation for the matching point, θ = θ_m, at which the slope of η_inner becomes equal to the constant slope of the η_outer, which prevails throughout most of the range 0 ≤ θ ≤ (π/2) one finds:

${{\eta }^{\prime }}_{\text{inner}}\left({\theta }_{\text{m}}\right)\approx -{{\eta }^{\prime }}_{\max }\Rightarrow {\theta }_{\text{m}}\approx \frac{\pi }{4}{\mu }^2\left(1+\frac{1}{2}{\mu }^2\right)$ . (58)

Requiring that η_inner and η_outer coincide at θ = θ_m yields:

${\eta }_{\text{outer}}\left(\theta \right)=\frac{2}{\pi }\left(\frac{\pi }{2}-\theta \right)+\frac{{\mu }^2}{4}\left(1-\frac{8}{{\pi }^2}\theta \right)+O\left({\mu }^4\right)$ . (59)

δ-function profile for acceleration As in previous sections, η''(θ) tends to a periodic sequence of δ-functions with alternating signs (Equation (18)). Again, this is well approximated by the cumulative change of η'(θ) over “the inner solution” as the turning point is crossed. Using Equations (55) and (58), one obtains the total change in the slope of the “outer solution” as the turning point is crossed:

$\begin{array}{c}{\eta }^{\prime }\left({\theta }_m\right)-{\eta }^{\prime }\left(-{\theta }_m\right)\approx {\displaystyle \underset{-{\theta }_m}{\overset{+{\theta }_m}{\int }}{{\eta }^{″}}_{inner}\left(\theta \right)\text{d}\theta }=-\frac{4\sqrt{2}}{\pi \mu }\sin \left(\mu \left(2+{\mu }^2\right)/\left(2\sqrt{2}\right)\right)\\ =-\frac{4}{\pi }\left(1+\frac{5}{12}{\mu }^2+O\left({\mu }^4\right)\right)\underset{\mu \to 0}{\to }-\frac{4}{\pi }\end{array}$ . (60)

In the extreme relativistic limit (μ → 0), the preceding analysis ensures that the linear “outer solution” describes the solution very well and the time interval around the turning points, where the solution is described by x_inner tends to zero. As a result, the asymptotic profile of the solution is also a saw-tooth one, as in Section 2.1.

2.3. Polynomial Nonlinearity: Coefficients That Grow with Energy

A new vista opens up when coefficients of the non-leading monomials in Equation (34) grow with energy, i.e., if for some p < N,

$\begin{array}{l}{a}_{2p+1}\underset{a\left(E\right)\to \infty }{\to }\alpha a{\left(E\right)}^{2\left(N-p\right)}\Rightarrow \\ a_{2p+1}{x}^{2p+1}/a_{2N+1}a{\left(E\right)}^{2N+1}\underset{x\to a\left(E\right)}{\to }{a}_{2p+1}{a}^{2p+1}/a_{2N+1}a{\left(E\right)}^{2N+1}\to O\left(1\right)\end{array}$ , (61)

if, in addition, a_2p+1 > 0 and a_2N+1 < 0. More than one term in the polynomial participates in nonlinear violence, and owing to the opposite signs, they compete with one another. For a → ∞, an additional leading monomial shows up in Equation (38):

${\eta }^{″}\left(\theta \right)+{\left(\frac{G}{2\pi }\right)}^2\left\{\alpha \eta {\left(\theta \right)}^{2p+1}-\eta {\left(\theta \right)}^{2N+1}\right\}=0$ . (62)

Unlike the case of constant coefficients in a polynomial potential (Section 2.1.3), the additional requirement of N → ∞ is not necessary here. The characteristics of the solutions for different values of p and N are similar, provided α = 1. Hence, one example of physical significance is discussed, the Duffing equation with an energy-dependent liner term:

$\ddot{x}\left(t\right)+\left(1+\frac{1}{2}a{\left(E\right)}^2\right)x-\frac{1}{2}x{\left(t\right)}^3$ . (63)

Equation (63) was derived by an elegant transformation of the equation of motion of the relativistic harmonic oscillator, Equation (41), using proper time (time as measured on the moving mass rather than time measured by an external observer) and scaled, dimensionless coordinates [46] [47] [51] [52] [53] . The linear and cubic terms are comparable near the turning points, x = ±a(E). Employing Equation (4), the scaled equation is:

${\eta }^{″}\left(\theta \right)+\frac{T^2}{4{\pi }^2}\eta \left(\theta \right)+\frac{T^2}{4{\pi }^2}\frac{1}{2}{a}^2\left(\eta \left(\theta \right)-\eta {\left(\theta \right)}^3\right)=0,\left(\eta \left(0\right)=1,{\eta }^{\prime }\left(0\right)=0\right)$ . (64)

The period T is given by [26] [27] [28] [29] :

$T=4E\left(i\frac{a}{2}\right)=\frac{8}{\sqrt{a^2+4}}E\left(\frac{a}{\sqrt{a^2+4}}\right)=8\frac{\log \left(2a\right)}{a}-8\frac{\log \left(2a\right)-1}{a^3}+\cdots$ , (65)

where properties of E(k) (see Equation (50)) are found in [54] . Using Equation (65), Equation (64) can be expressed as a traditional boundary-layer problem:

$\begin{array}{l}0=\frac{1}{{\left(aT\right)}^2}{\eta }^{″}\left(\theta \right)+\frac{1}{4{\pi }^2{a}^2}\eta \left(\theta \right)+\frac{1}{8{\pi }^2}\left(\eta \left(\theta \right)-\eta {\left(\theta \right)}^3\right)\\ \approx \frac{1}{{\left(8\log \left(2a\right)\right)}^2}{\eta }^{″}\left(\theta \right)+\frac{1}{4{\pi }^2{a}^2}\eta \left(\theta \right)+\frac{1}{8{\pi }^2}\left(\eta \left(\theta \right)-\eta {\left(\theta \right)}^3\right),\eta \left(2\pi +\theta \right)=\eta \left(\theta \right)\end{array}$ (66)

The highest derivative is multiplied by a small parameter, which varies slowly as a → ∞.

Equation (66) has a closed-form solution [40] [41] [42] [43] [44] [47] , given by:

$\eta \left(\theta \right)=sn\left[\frac{\sqrt{a^2+4}}{2}\frac{T\left(\theta +\pi /2\right)}{2\pi },\frac{a}{\sqrt{a^2+4}}\right]$ . (67)

In Equation (67), sn is a Jacobi elliptic function [54] .

Applying energy conservation to Equation (66), or computing the derivative of η'(θ = π/2) through Equation (67), one obtains the maximal value of the derivative to be:

${{\eta }^{\prime }}_{\max }=\frac{\sqrt{4+a^2}}{4\pi }T=\frac{2}{\pi }\log \left(2a\right)+\frac{2}{a^2}\left(1+\log \left(2a\right)\right)$ . (68)

The “inner solution” describes the rapid change around θ = π/2. It is obtained by solving for η(θ) in a power series in (π/2 − θ): Through O((π/2 − θ)²), it is given by:

${\eta }_{\text{inner}}\left(\theta \right)={{\eta }^{\prime }}_{\max }\left(\frac{\pi }{2}-\theta \right)+A{\left(\frac{\pi }{2}-\theta \right)}^2$ . (69)

The “outer solution”, describing the moderate change around θ = 0 is found by solving Equation (66) around η(θ) = 1:

${\eta }_{\text{outer}}\left(\theta \right)=1-\delta \left(\theta \right)$ (70)

Substituting Equation (70) in Equation (66), the linearized equation for $\delta \left(\theta \right)\ll \text{1}$ becomes:

$\frac{T^2}{4{\pi }^2}\left\{1-\left(a^2-1\right)\delta \left(\theta \right)\right\}+{\delta }^{″}\left(\theta \right)=0,\delta \left(0\right)={\delta }^{\prime }\left(0\right)=0$ . (71)

Equation (71) is solved by

$\delta \left(\theta \right)=2\sinh {\left(\frac{\sqrt{a^2-1}T\theta }{4\pi }\right)}^2/\left(a^2-1\right)\approx \frac{2E\left(i\frac{A}{2}\right)}{{\pi }^2}{\theta }^2$ . (72)

The matching point, θ_m, is found by requiring

${\eta }_{\text{inner}}\left({\theta }_{\text{m}}\right)={\eta }_{\text{outer}}\left({\theta }_{\text{m}}\right),{{\eta }^{\prime }}_{\text{inner}}\left({\theta }_{\text{m}}\right)={{\eta }^{\prime }}_{\text{outer}}\left({\theta }_{\text{m}}\right)$ . (73)

One finds:

$\begin{array}{l}{\theta }_{\text{m}}=\frac{\pi }{2}-\frac{\pi }{2\log \left(2a\right)}+2\pi \frac{\log {\left(2a\right)}^2-\log \left(2a\right)+1}{a^2}+\cdots \\ A={{\eta }^{\prime }}_{\max }^2\left(\frac{1}{4}+\frac{\log {\left(2a\right)}^2-4\log \left(2a\right)}{2a^2}\right)\end{array}$ . (74)

The rate of approach as a → ∞ is slow. The matching point, θ_m, deviates from (π/2) by O(1/log2a). Figure 5 shows the explicit solution for η (Equation (67)), η_inner Equations (69) and (74), η_outer Equations (70) and (72) and the location of θ_m, Equation (74).

3. V(x) Not Symmetric under Reflection in x, (1 + 1) Dimensions

3.1. General Comments

If V(x) is not symmetric under reflection in x, the amplitudes at the two turning points are not equal up to a sign. Consequently, the analysis has to be performed separately close to any turning point, around which nonlinear violence prevails. The overall scaling of the solution (x = a ∙ η) in Equation (4) has to be replaced near each turning point by (x = x_tur(E)×η), where x_tur(E) is the amplitude at that turning point. Finally, while the transformation (t → θ) of Equation (4) yields

Figure 5. Solution of Equation (64). Full line: closed-form solution (Equation (67)); dotted line: outer solution (Equations (70) and (74)); dashed line: inner solution (Equations (69) and (72)).

(2π)-periodicity in θ, there is no symmetry between (0 ≤ θ ≤ π) and (π ≤ θ ≤ 2π). For example, if only one turning point exhibits nonlinear violence, T is a sum of two terms: a vanishingly small interval of time (spent around that turning point affected by nonlinear violence) and the rest of the time, spent over most of the oscillation.

Except for technical details, the analysis, specifically, the procedure for matching the “inner” and “outer” solutions, is the same as presented in Sections 2.1 - 2.4. Hence, instead of analyzing a general system, the analysis is carried out in Sections 3.2 - 3.3 for a specific example, where the solution is affected by nonlinear violence only around one turning point.

3.2. The System

Consider a harmonic oscillator with an exponential nonlinearity:

$\ddot{x}\left(t\right)+x\left(t\right)+{\text{e}}^{x\left(t\right)}=0$ . (75)

The potential energy is asymmetric:

$\begin{array}{l}V\left(x\right)=\frac{1}{2}{x}^2+{\text{e}}^x\Rightarrow \\ V\left(x\right)\underset{x\to -\infty }{\to }\frac{1}{2}{x}^2,V\left(x\right)\underset{x\to +\infty }{\to }{\text{e}}^x\\ E=\frac{1}{2}{x}_{\max }{\left(E\right)}^2+{\text{e}}^{x_{\max }\left(E\right)}\end{array}$ . (76)

For large E, at large negative x, V(x) of Equation (76) tends to that of a harmonic oscillator. Hence, the profile of oscillations is expected to be close to a trigonometric function. Nonlinear violence occurs only around the turning point, x_max(E) > 0, as there the potential energy grows rapidly. (x_max(E) → ∞ as E → ∞.) If the initial conditions are such the x(0) is that turning point, so that $\stackrel{\dot{}}{x}\left(0\right)$ = 0, turning points occur at t₀ = 0, T, 2T, …

As E → ∞, for x(t) < 0, the exponential term in Equation (75) becomes rapidly negligible. Hence, the time spent over the x(t) < 0 branches of oscillations approaches π (half of a full sinusoidal oscillation). The time spent around the turning point, which is governed by nonlinear violence, is vanishingly small. Hence, T tends to π as E → ∞. For x(0) = 10, numerical integration yields T = 3.25.

Numerical Examples

The boundary-layer structure that emerges around x_max(E) is demonstrated Figure 6 and Figure 7 for an initial amplitude, x(0) = x_max(E) = 10. The details around a turning point are shown in Figure 8 and Figure 9. Near the turning point, x(t) is highly curved, soon evolving into a profile that is well approximated by linear time dependence. Over the same region, $\stackrel{\dot{}}{x}\left(t\right)$ varies rapidly between two constants, ±( $\sqrt{2E}$ ). This boundary-layer behavior is soon replaced by a trigonometric profile that dominates most of the period of oscillations.

Figure 6. x(t), solution of Equation (75). x(0) = 10, $\stackrel{\dot{}}{x}\left(0\right)=0$ .

Figure 7. $\stackrel{\dot{}}{x}\left(t\right)$ , solution of Equation (75). x(0) = 10, $\stackrel{\dot{}}{x}\left(0\right)=0$ .

Figure 8. x(t), solution of Equation (75) near positive turning point. x(0) = 10, $\stackrel{\dot{}}{x}\left(0\right)=0$ .

Figure 9. $\stackrel{\dot{}}{x}\left(t\right)$ , solution of Equation (75) near positive turning point. x(0) = 10, $\stackrel{\dot{}}{x}\left(0\right)=0$ .

3.3. Nonlinear Violence Analysis

To find an approximation for the solution of Equation (75) near the turning point, a = x_max(E) > 0, for large E, write x(t) as:

$x\left(t\right)=a\eta \left(\tau \right),\tau =\alpha t,0\le \eta \le 1$ . (77)

Equation (75) becomes

$\frac{1}{a}{\text{e}}^{a\eta \left(\tau \right)}+\eta \left(\tau \right)+{\alpha }^2{\eta }^{″}\left(\tau \right)=0$ . (78)

Energy conservation for Equation (75) yields

$a\underset{E\gg 1}{=}\log E-\frac{{\left(\log E\right)}^2}{2E}+o\left(\frac{{\left(\log E\right)}^2}{2E}\right)$ . (79)

To find α, substitute Equation (77) in Equation (75) and approximate the “inner solution” by:

${\eta }_{\text{inner}}\left(\tau \right)=1-{\tau }^2+O\left({\tau }^4\right)$ . (80)

To leading order, one obtains

$\alpha =\sqrt{\frac{E}{2\log E}}\left(1-\frac{{\left(\log E\right)}^2}{4E}+\frac{3\log E}{4E}+o\left(\frac{1}{E}\right)\right)$ . (81)

As the solution evolves into a trigonometric profile once away from the turning point, we expect the “outer solution” to provide a good approximation (linear in time) for the full solution for a limited range in time. For this reason, scaling of the time variable with α of Equation (81), and a of Equation (79) will be employed also in the generation of the “outer solution”. Equation (75) becomes:

$2\left(\frac{\log E}{E}+E^{\eta \left(\tau \right)-1}\right)\left(1+O\left(1/E\right)\right)\eta \left(\tau \right)+{\eta }^{″}\left(\tau \right)=0$ . (82)

For E → ∞ and (1 − |η(τ)|) = O(1)) > 0, Equation (82) tends to:

${\eta }^{″}\left(\tau \right)=0$ . (83)

This yields the “outer solution”:

${\eta }_{\text{outer}}\left(\tau \right)=A+B\tau$ . (84)

Matching of the “inner” and “outer” solutions by making their values and their derivatives equal at t = t_m yields the matching point:

$t_{\text{m}}=\sqrt{\frac{2}{E}}\left(1+\frac{\log E}{2E}+o\left(\frac{1}{E}\right)\right)$ . (85)

For a turning point at t = t₀, the “inner solution” is found to be

$x{\left(t\right)}_{\text{inner}}\approx a\left(1-\frac{E{\left(t-t_0\right)}^2}{2\log E}+\frac{E^2{\left(t-t_0\right)}^4}{24\log E}\right)+O\left(E^3{\left(t-t_0\right)}^6\right)$ , (86)

and the outer solution:

$x{\left(t\right)}_{\text{outer}}\approx a\left(\left(1+\frac{7}{6\log E}+O\left(\frac{1}{E}\right)\right)+\frac{\sqrt{2E}}{\log E}\left(1+\frac{\log E}{2E}-O\left(\frac{1}{E}\right)\right)\left|t-t_0\right|\right)$ . (87)

A comparison of the approximations provided by the “inner solution” and the “outer solution” with the numerical solution of Equation (75) around the turning point amplitude, x(t₀) = x_max = 10 and t₀ = 3.25, is shown in Figure 10. The “inner solution” well approximates the curved part of the full solution near the turning point and the “outer solution” does so (for a while) away from the turning point.

Velocity change across turning point Using Equations (83) and (84), the cumulative change of η'(t) over the diminishingly small interval, (t₀ − t_m, t₀ + θ_m,), (θ_m → 0 for E → ∞) is found to coincide with the change in the slope of the “outer solution”, x_out, as one crosses the turning point from left to right:

Figure 10. Comparison of quadratic and quartic approximations for “inner solution” (Equation (86), dashed and dotted lines, respectively) and “outer solution” (Equation (87), dot-dashed) with numerical solution of Equation (75) near t₀ = 3.25 turning point. x(0) = 10, $\stackrel{\dot{}}{x}\left(0\right)=0$ .

$x^{\prime }\left(t_0+t_m\right)-x^{\prime }\left(t_0-t_m\right)={\displaystyle \underset{t_0-t_m}{\overset{t_0+t_m}{\int }}{x^{″}}_{inner}\left(t\right)\text{d}t}=-2\frac{\sqrt{2E}}{\log E}\left(1+O\left(1/E\right)\right)$ , (88)

again, asymptotically yielding the overall change in x'(t) as the turning point is crossed.

4. (1 + 2) Dimensions: Motion in Central Force with Non-Zero Angular Momentum

At high E, the motion of a particle in a central field in two space dimensions offers new interesting boundary-layer phenomena if the angular momentum is not equal to zero. The evolution of the radius, r, of a unit mass particle is determined by [36] :

$\begin{array}{l}{r}^{″}+{\partial }_r{V}_{eff}\left(r\right)=0\\ V_{eff}\left(r\right)=l^2/\left(2r^2\right)+V\left(r\right)\end{array}$ . (89)

Here V(r) is the potential and l is the conserved angular momentum:

$l=r^2\left(t\right)\stackrel{\dot{}}{\phi }\left(t\right)=Const.$ (90)

φ is the angle of rotation of the moving mass around the origin.

4.1. Gravitational Potential

4.1.1. Analysis of Solution

For E < 0, the solution classical motion under a gravitational potential,

$V\left(r\right)=-1/r$ , (91)

is bounded and periodic. For E > 0, the solution visits the point of the minimal radius, r_min(E), once and wonders off to infinity [36] . Boundary-layer characteristics emerge in the vicinity of the smallest radius, r_min, when $E\gg \text{1}$ if l ≠ 0.

4.1.2. Turning Point at r_min

Energy conservation for Equation (89) with V(r) of Equation (91),

$E=l^2/\left(2r{\left(t\right)}^2\right)-1/r\left(t\right)+\frac{1}{2}\stackrel{\dot{}}{r}{\left(t\right)}^2=l^2/\left(2r_{\text{min}}{\left(E\right)}^2\right)-1/r_{\text{min}}\left(E\right)$ , (92)

yields:

$r_{\text{min}}\left(E\right)=\frac{\sqrt{2l^{2}E+1}-1}{2E}=\frac{l}{\sqrt{2E}}-\frac{1}{2E}+O\left(E^{-3/2}\right)$ (93)

and

$\stackrel{\dot{}}{r}\left(t\right)=\pm \sqrt{2E+\left(2/r\left(t\right)\right)-{\left(l/r\left(t\right)\right)}^2}\underset{r\left(t\right)\to \infty }{\to }\pm \sqrt{2E}$ (94)

Hence, the “outer solution” is given by:

$r_{\text{outer}}\left(t\right)=A+\sqrt{2E}t$ (95)

For the “inner solution”, as r_min is a local minimum around t = 0, write:

$r_{\text{inner}}\left(t\right)=r_{\min }\left(1+\alpha t^2\right)$ . (96)

Substituting Equation (96) in Equation (89) yields:

$\alpha =\frac{l^2}{2r_{\min }{\left(E\right)}^4}\left(1-\left(r_{\min }\left(E\right)/l^2\right)+O\left(r_{\min }{\left(E\right)}^2\right)\right)$ . (97)

Using Equations (95) - (97), matching the “inner” and “outer” solutions and their time derivatives at the matching time, t_m, yields:

$\begin{array}{l}{t}_m=\frac{\sqrt{2E}{r}_{\text{min}}{\left(E\right)}^3}{l^2-r_{\text{min}}\left(E\right)}=\frac{l}{2E}\left(1+O\left(\frac{1}{\sqrt{E}}\right)\right),\\ A=\frac{l^2}{2\sqrt{2E{l}^2+1}}=\frac{l}{2\sqrt{2E}}\left(1+O\left(1/E\right)\right).\end{array}$ (98)

Thus, the matching point, t_m, vanishes like E^-1.

Velocity change across turning point, r_min Using Equations (90) and (93) - (98), the total change in the velocity $\stackrel{\dot{}}{r}\left(t\right)$ , around r_min over the time interval −t_m ≤ t ≤ + t_m, computed in terms of r_inner, reproduces the full change in the slope of the “outer solution” from negative to positive linear slopes:

$\stackrel{\dot{}}{r}\left(+t_m\right)-\stackrel{\dot{}}{r}\left(-t_m\right)\approx {\displaystyle \underset{-t_m}{\overset{+t_m}{\int }}{\ddot{r}}_{\text{inner}}\left(t\right)\text{d}t}=2r_{\text{min}}\left(E\right)\alpha t_m=2\sqrt{2E}\left(1+O\left(\frac{1}{\sqrt{E}}\right)\right)$ . (99)

Rapid rotation at turning point, r_min

The last interesting aspect is the variation of the angle of rotation, φ, in time as E → ∞. Close to the turning point, Equations (90), (94), (96) and (97) yield:

$\begin{array}{l}\stackrel{\dot{}}{\phi }\left(t\to 0\right)=\frac{l}{r_{\text{min}}{\left(E\right)}^2}=\frac{2E}{l}\left(1+\left(\sqrt{2/E}\right)/l+O\left(1/E\right)\right)\\ \stackrel{\dot{}}{\phi }\left(t\to t_m\right)=\frac{l}{r{\left(t_m\right)}^2}=\frac{8E}{9l}\left(1+\left(4/3\right)\left(\sqrt{2/E}\right)/l+O\left(1/E\right)\right)\end{array}$ . (100)

Thus, $\stackrel{\dot{}}{\phi }\left(t\right)$ becomes infinite in the vicinity of the turning point as E → ∞. However, soon the “outer solution” takes over, for which:

$\stackrel{\dot{}}{\phi }\left(t\gg t_m\right)=\frac{l}{r_{\text{outer}}{\left(t\right)}^2}\approx \frac{1}{2E{t}^2}\underset{E\to \infty }{\to }0$ . (101)

In summary, close to r_min, the angle, φ, varies rapidly between two constants over a vanishingly small time interval around t = 0. Away from that interval it varies very slowly with time up to the large turning point, r_max. In the limit E → ∞, φ exhibits a finite jump around the r_min. This point will be expanded in Section (4.2).

4.2. Positive Potential That Grows as r → ∞

4.2.1. Analysis of Solution

Consider in Equation (89) V(r), which grows rapidly at large r. To be specific choose:

$V\left(r\right)={\text{e}}^r$ . (102)

The motion is oscillatory as the total energy is:

$E=\left(l^2/\left(2r^2\right)\right)+{\text{e}}^r$ . (103)

Nonlinear violence, leading to a boundary-layer phenomenon, emerges around the turning point at the maximal radius, r_max. For $E\gg \text{1}$ :

$r_{\text{max}}\approx \log E$ . (104)

Nonlinear violence emerges also in the vicinity of r_min if l ≠ 0 and the contribution of the angular momentum term dominates:

$\frac{l^2/\left(2r_{\min }^2\right)}{E}\underset{E\to \infty }{\to }1\Rightarrow \frac{V\left(r_{\min }\right)}{l^2/\left(2r_{\min }^2\right)}\underset{E\to \infty }{\to }0$ . (105)

Equation (105) holds for V(r) of Equation (102). The analysis is identical to that presented in Section 4.1, hence, not repeated here.

Numerical examples are shown for E = 10⁶ and l = 1. Figure 11 shows the numerical solution for r(t) of Equation (89) with the potential of Equation (102). The profile is almost linear in t, except for rapid bending close to the turning points. Figure 12 shows the numerical solution for $\stackrel{\dot{}}{r}\left(t\right)$ . The profile is almost constant, except for rapid change of sign close to the turning points. The rapid transition of r(t) from the highly curved “inner solution” to the linear profile of the “outer solution” near r_min and r_max is demonstrated in Figure 13 and Figure 14, respectively.

4.2.2. Turning Point at r_max

Energy conservation for Equation (89) with V(r) of Equation (102),

$E=l^2/\left(2r{\left(t\right)}^2\right)+{\text{e}}^{r\left(t\right)}+\frac{1}{2}\stackrel{\dot{}}{r}{\left(t\right)}^2=l^2/\left(2r_{\text{max}}{\left(E\right)}^2\right)+{\text{e}}^{r_{\max }}$ , (106)

yielding Equation (104) for r_max and

Figure 11. Numerical solution of Equation (89) for r(t) with exponential potential energy (Equation (102)), E = 10⁶. Dashed lines represent r_min and r_max.

Figure 12. Numerical solution of Equation (89) for $\stackrel{\dot{}}{r}\left(t\right)$ with exponential potential energy (Equation (102)), E = 10⁶ and l = 1. Dashed lines represent the constant values of $\stackrel{\dot{}}{r}\left(t\right)$ away from turning points.

Figure 13. Numerical solution of Equation (89) for r(t) with exponential potential energy (Equation (102)) in vicinity of r_min, E = 10⁶ and l = 1.

Figure 14. Numerical solution of Equation (89) for r(t) with exponential potential energy (Equation (102)) in vicinity of r_max, E = 10⁶ and l = 1.

$\stackrel{\dot{}}{r}\left(t\right)=\pm \sqrt{2E-2{\text{e}}^{r\left(t\right)}-{\left(l/r\left(t\right)\right)}^2}\underset{\begin{array}{l}E\to \infty \\ r\left(t\right)\ll r_{\max }\end{array}}{\to }\pm \sqrt{2E}$ . (107)

Hence, the “outer solution” is given by:

$r_{\text{outer}}\left(t\right)=A+\sqrt{2E}t$ . (108)

For the “inner solution”, as r_max is a local minimum around t = 0, write:

$r_{\text{inner}}\left(t\right)=r_{\max }\left(1-\alpha t^2\right)$ . (109)

Substituting Equation (109) in Equation (89) with the potential of Equation (102) yields:

$\alpha =\frac{{\text{e}}^{r_{\max }}}{2r_{\max }}\left(1-l^2{\text{e}}^{-r_{\max }}/r_{\max }^3\right)\approx \frac{E}{2\log E}\left(1-l^2/\left(E{\left(\log E\right)}^3\right)\right)$ . (110)

Using Equations (108)-(110), matching the “inner” and “outer” solutions and their time derivatives at the matching time, t_m, yields:

$\begin{array}{l}{t}_m=\frac{\sqrt{2}{\text{e}}^{-r_{\max }/2}}{1-\left({\text{e}}^{-r_{\max }}{l}^2/r_{\max }^3\right)}=\frac{\sqrt{2}}{E^{/12}}\left(1+O\left(\frac{1}{E}\right)\right)\\ A=\left(1+r_{\max }\right)\left(1+O\left(\frac{{\text{e}}^{-r_{\max }}{l}^2}{r_{\max }^3}\right)\right)\approx 1+\log E\left(1+O\left(\frac{l^2}{E{\left(\log E\right)}^3}\right)\right)\end{array}$ . (111)

Thus, the matching point, t_m, vanishes as E^-1/2.

Velocity change across turning point, r_max Using Equations (108) and (111), the total change in the velocity $\stackrel{\dot{}}{r}\left(t\right)$ around r_max as one crosses the turning point, is well reproduced by the contribution of r_inner:

$\stackrel{\dot{}}{r}\left(+t_m\right)-\stackrel{\dot{}}{r}\left(-t_m\right)={\displaystyle \underset{-t_m}{\overset{+t_m}{\int }}{\ddot{r}}_{\text{inner}}\left(t\right)\text{d}t}=-2\sqrt{2{\text{e}}^{r_{\max }}}\approx -2\sqrt{2E}$ . (112)

Again, this is the total change in ${\stackrel{\dot{}}{r}}_{\text{outer}}\left(t\right)$ .

4.2.3. Variation of Rotation Angle near Turning Points

Near r_min

The angular velocity, $\stackrel{\dot{}}{\phi }\left(t\right)$ , is shown in Figure 15 for the numerical solution of Equation (89) with the potential of (102). It is highly peeked near r_min and negligibly small elsewhere. Details are shown in Figure 16 and Figure 17 where $\stackrel{\dot{}}{\phi }\left(t\right)$ is plotted in the vicinity of each turning point. Using Equations (93) and (96)-(98), the change of $\phi \left(t\right)$ when crossing r_min over the time interval (−t_m, +t_m) is found to be:

Figure 15. Computation of $\stackrel{\dot{}}{\phi }\left(t\right)$ using numerical solution of Equation (89) with exponential potential energy (Equation (102)), E = 10⁶ and l = 1.

Figure 16. Computation of $\stackrel{\dot{}}{\phi }\left(t\right)$ using numerical solution of Equation (89) with exponential potential energy (Equation (102)) in vicinity of r_min, E = 10⁶ and l = 1.

Figure 17. Computation of $\stackrel{\dot{}}{\phi }\left(t\right)$ using numerical solution of Equation (89) with exponential potential energy (Equation (102)) in vicinity of r_max, E = 10⁶ and l = 1.

$\Delta \phi ={\displaystyle {\int }_{-t_{\text{m}}}^{+t_{\text{m}}}\left(L/r_{\text{inner}}{\left(t\right)}^2\right)\text{d}t}=\frac{2}{3}+2\tan \left(1/\sqrt{2}\right)=1.53709$ (113)

If instead of r_inner we use the numerical solution of Equation (89) for r(t), with l = 1, the potential of Equation (102), and E = 10⁶ in Equation (113) the computed jump obtained is 1.5715. The exact result is that, as E → ∞, Δφ jumps across r_min by (π/2) = 1.5708 [36] . Both results are very good approximations deviating from (π/2) owing to the approximations employed in expressing r_inner, r_outer and t_m.

Near r_max

Assuming that r_max occurs at t = 0, one finds in the vicinity of r_max:

$\begin{array}{l}{r}_{\text{inner}}\left(t\right)=r_{\text{max}}\left(1-\left({\text{e}}^{r_{\text{max}}}/\left(2r_{\text{max}}\right)\right)\left(1+O\left(\left({\text{e}}^{-r_{\text{max}}}/r_{\max }^3\right)\right)\right)t^2\right)\\ r_{\text{outer}}=A+\sqrt{2E}t,A=r_{\max }+1+O\left({\text{e}}^{-r_{\text{max}}}\left(1/r_{\max }^2\right)\right)\\ t_{\text{m}}=\sqrt{2}{\text{e}}^{-r_{\max }/2}\left(1+O\left(\left({\text{e}}^{-r_{\text{max}}}/r_{\max }^2\right)\right)\right)\end{array}$ . (114)

The integral of Equation (113) for r_inner of Equation (114) reveals that the change in φ(t) across r_max in the time interval (−t_m, + t_m), with t_m of Equation (114) vanishes as E → ∞:

$\Delta \phi =O\left({\text{e}}^{-r_{\max }/2}/r_{\max }^2\right)=O\left(1/\left(E^{1/2}{\left(\log E\right)}^2\right)\right)$ , (115)

In summary, in agreement with what the numerical solution shows, φ(t) varies rapidly (a jump of (π/2) in the limit E → ∞) when crossing r_min and by a vanishingly small change when crossing r_max. Outside the short time intervals spent around r_min (Equation (93)) and r_max (Equation (104)), the total variation of $\phi \left(t\right)$ along each branch of the oscillation cycle is (3/4π), occurring at a very low rate, O(1/(E ∙ t²)).

Comment: Choice of the exponential potential of Equation (102), the contribution of which to V_eff of Equation (89) becomes negligible at small r, was motivated by the role played by the angular momentum in generating rapid rotation around r_min. A detailed analysis (not presented here) shows that, if the exponential potential of Equation (89) is replaced by V(r) such that its contribution to V_eff of Equation (89) dominates also at small r, then, depending on the details of V(r), rapid rotation around r_min may or may not emerge.

5. Concluding Comments

Nonlinear oscillations with high oscillation energy are of relevance in many areas of science and technology: communications science, electrical circuits, electronic devices, biological systems and more. A recent example has emerged in low-temperature physics, in the oscillations of ions in a Bose-Einstein condensate [55] , where the speed of light is smaller than in vacuum by orders of magnitude and a saw-tooth profile for the oscillations has been observed.

There is a large family of oscillatory systems whose solutions are affected by nonlinear violence, leading to boundary layer-like characteristics when the oscillation energy is extremely high. These characteristics emerge independently of whether the dynamical equation of the oscillator can or cannot be cast in the traditional form of boundary-layer problems.

Planar motion of a particle in a central field displays a new unique phenomenon, which opens the door to interesting behavior of oscillators in higher space dimensions. Boundary-layer characteristics emerge around the turning point of the largest radius owing to the growth of the potential energy. If the angular momentum, l, does not vanish, nonlinear violence may emerge around the turning point of the smallest radius, r_min. The latter is accompanied by a boundary-layer characteristic behavior of the angle that determines the position of the rotating particle in the plane. At large overall energy, the change in the angle when crossing r_min is very rapid and tends to a discrete jump of (π/2) as E → ∞.

In all the cases discussed here, the change of velocity and of rotation angle when a turning point is crossed are well approximated by computing them with the relevant “inner solution” over the time interval: (−t_m ≤ t ≤ +t_m).

In the case of (1 + 1)-dimensional energy conserving nonlinear oscillatory systems, the effect of nonlinear violence depends on the detailed structure of the equations of motion. If the potential energy of the system is not symmetric under space reflection, nonlinear violence may emerge around one or both turning points. If the potential energy is symmetric under reflection in x, these characteristics emerge at both turning points. The solution for x(t) then tends to a simple periodic profile: examples of saw-tooth and step-function profile were discussed.

The analysis discussed in this paper leads to the conclusion that, as E → ∞, the asymptotic profiles of solutions of a large family of oscillatory systems tend to similar standard forms, such as saw-tooth or step-function profiles. Hence, if such profiles are observed experimentally in oscillating physical systems, it may be difficult to decipher the equations of motion. A possible future research direction that will improve our understanding of high-energy nonlinear oscillations is the development of an expansion analysis around the asymptotic (saw-tooth or step-function) profiles, for the deviations of the actual solution from the asymptotic profile. This may be very important in deciphering the equations of motion governing the dynamical system, for which only experimental observations exist and theory has not been developed yet.

Author Responsibility

1) The author made all contributions to the conception or design of the work. No data acquisition was required. No new software was created;

2) The author drafted the work or revised it critically for important intellectual content;

3) The author approved the version to be published;

4) The author agrees to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Conflicts of Interest

The author declares that he has no conflict of interest.

References

[1]	Poincaré, H. (1993) New Methods of Celestial Mechanics. Springer-Verlag, Berlin.
[2]	Lindstedt, A. (1884) Ueber die Bestimmung der gegenseitigen Entfernungen in dem Probleme der drei Körper. Astronomische Nachrichten, 107, 197-214. https://doi.org/10.1002/asna.18841071301
[3]	Guckenheimer, J. and Holmes, P. (1983) Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer, Berlin. https://doi.org/10.1007/978-1-4612-1140-2
[4]	Arnold, V.I. (1989) Mathematical Methods of Classical Mechanics. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4757-2063-1
[5]	Kahn, P.B. and Zarmi, Y. (1998) Nonlinear Dynamics: Exploration through Normal Forms. Wiley, Hoboken.
[6]	Giampaolo, C. and Gaeta, G. (1999) Symmetry and Perturbation Theory in Nonlinear Dynamics. Springer, Berlin.
[7]	Murdock, J. (2003) Normal Forms and Unfoldings If Local Dynamical Systems. Springer, Berlin. https://doi.org/10.1007/b97515
[8]	Chow, S.N. (2010) Normal Forms and Bifurcations of Planar Vector Fields. Cambridge University Press, Cambridge.
[9]	Nayfeh, A.H. (2011) The Method of Normal Forms. Wiley, Hoboken. https://doi.org/10.1002/9783527635801
[10]	Cabral, H.E. and Dias, L.B. (2023) Normal Forms and Stability of Hamiltonian Systems. Springer, Berlin. https://doi.org/10.1007/978-3-031-33046-9
[11]	Mitropolsky, Iu.A. (1967) Averaging Method in Nonlinear Mechanics. International Journal of Non-Linear Mechanics, 2, 69-96. https://doi.org/10.1016/0020-7462(67)90020-0
[12]	Sanders, J. and Verhulst, F. (1985) Averaging Methods in Nonlinear Dynamical Systems. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4757-4575-7
[13]	Murdock, J.A. (1991) Perturbations: Theory and Methods. Wiley, Hoboken.
[14]	Nayfeh, A.H. (1985) Problems in Perturbation. Wiley, Hoboken.
[15]	Nayfeh, A.H. (1993) Method of Normal Forms. Wiley, Hoboken.
[16]	Nayfeh, A.H. (1993) Introduction to Perturbation Technique. John Wiley & Sons, Hoboken.
[17]	Nayfeh, A.H. and Mook, D.T. (1995) Nonlinear Oscillations. Wiley, Hoboken. https://doi.org/10.1002/9783527617586
[18]	Nayfeh, A.H. (2000) Perturbation Methods. Wiley, Hoboken. https://doi.org/10.1002/9783527617609
[19]	He, J.H. (2000) Α Review on Some New Recently Developed Nonlinear Analytical Techniques. The International Journal of Nonlinear Sciences and Numerical Simulation, 1, 51-70. https://doi.org/10.1515/IJNSNS.2000.1.1.51
[20]	Cveticanin, L. (2014) Strongly Nonlinear Oscillators. Springer, Berlin. https://doi.org/10.1007/978-3-319-05272-4
[21]	Ismail, G.M. (2021) Highly Accurate Analytical Solution for Free Vibrations of Strongly Nonlinear Duffing Oscillator. Journal of Low Frequency Noise, Vibration and Active Control, 41, 223-229. https://doi.org/10.1177/14613484211034009
[22]	Tian, Y. and Feng, Q.Q. (2022) A Short Review on Approximate Analytical Methods for Nonlinear Problems. Thermal Science, 26, 2607-2628. https://doi.org/10.2298/TSCI2203607T
[23]	Beléndez, A., et al. (2008) Higher-Order Approximate Solutions to the Relativistic and Duffing-Harmonic Oscillators by Modified He’s Homotopy Methods. Physica Scripta, 77, Article ID: 025004. https://doi.org/10.1088/0031-8949/77/02/025004
[24]	Rosenberg, R.M. (1963) The Ateb(h)-Functions and Their Properties. Quarterly of Applied Mathematics, 21, 34-37. https://doi.org/10.1090/qam/143948
[25]	Mishchenko, E. and Rozov, N.Kh. (1980) Differential Equations with Small Parameters and Relaxation Oscillations. Plenum Press, New York. https://doi.org/10.1007/978-1-4615-9047-7
[26]	Pilipchuk, V.N. (1985) The Calculation of Strongly Non-Linear Systems Close to Vibration Impact Systems. Journal of Applied Mathematics and Mechanics, 49, 572-578. https://doi.org/10.1016/0021-8928(85)90073-5
[27]	Pilipchuk, V.N. (2010) Nonlinear Dynamics: Between Linear and Impact Limits. Springer, New York.
[28]	Gendelman, O., et al. (1999) Transition from Strongly to Weakly Nonlinear Motions of Damped Nonlinear Oscillators. Nonlinear Dynamics, 20, 99-114. https://doi.org/10.1023/A:1008354208466
[29]	Gendelman, O. and Vakakis, A.F. (2000) Transitions from Localization to Non-Localization in Strongly Nonlinear Damped Oscillators. Chaos, Solitons and Fractals, 11, 1535-1542. https://doi.org/10.1016/S0960-0779(99)00076-4
[30]	Awrejcewicz, J. and Andrianov, I.V. (2002) Oscillations of Non-Linear System with Restoring Force Close to Sign (x). Journal of Sound and Vibration, 252, 962-966. https://doi.org/10.1006/jsvi.2001.3666
[31]	Cveticanin, L. (2009) The Approximate Solving Methods for the Cubic Duffing Equation Based on the Jacobi Elliptic Functions. International Journal of Nonlinear Sciences and Numerical Simulation, 10, 1491-1516. https://doi.org/10.1515/IJNSNS.2009.10.11-12.1491
[32]	Andrianov, I.V., et al. (2011) Periodical Solutions of Certain Strongly Nonlinear Wave Equations. Numerical Analysis and Applicable Mathematics, 1389, 442-445. https://doi.org/10.1063/1.3636758
[33]	Andrianov, I.V., et al. (2014) Asymptotic Methods in the Theory of Plates with Mixed Boundary Conditions. Wiley, Chichester. https://doi.org/10.1002/9781118725184
[34]	Andrianov, I.V., et al. (2016) Analytical Approximation of Periodic Ateb-Functions via Elementary Functions. AIP Conference Proceedings, 1773, Article ID: 040001. https://doi.org/10.1063/1.4964964
[35]	Andrianov, I.V., et al. (2016) Analytical Approximation of Periodic Ateb-Functions via Elementary Functions. Proceedings 5th International Conference on Nonlinear Dynamics, Kharkov, 27-30 September 2016.
[36]	Goldstein, H. (1981) Classical Mechanics. Addison-Wesley, Boston, Ch. 3.
[37]	Guckenheimer, J. and Holmes, P. (1983) Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4612-1140-2
[38]	Mickens, R.E. (1981) Nonlinear Oscillations. Cambridge University Press, New York.
[39]	Rand, R.H. (2012) Lecture Notes on Nonlinear Vibrations. Cornell University, Ithaca, 13-17.
[40]	Salas, A.H. and Castillo, J.E. (2014) Exact Solutions to Cubic Duffing Equation for a Nonlinear Electrical Circuit. Visión Tecnológica, 7, 46-53.
[41]	Maimistov, A.I. (2003) Propagation of an Ultimately Short Electromagnetic Pulse in a Nonlinear Medium Described by the Fifth-Order Duffing Model. Optics and Spectroscopy, 94, 251-257. https://doi.org/10.1134/1.1555186
[42]	Elías-Zúniga, A. (2013) Exact Solution of the Cubic-Quintic Duffing Oscillator. Applied Mathematical Modelling, 37, 2574-2579. https://doi.org/10.1016/j.apm.2012.04.005
[43]	Beléndez, A., et al. (2016) Exact Solution for the Unforced Duffing Oscillator with Cubic and Quintic Nonlinearities. Nonlinear Dynamics, 86, 1687-1700. https://doi.org/10.1007/s11071-016-2986-8
[44]	Beléndez, A., et al. (2017) Closed-Form Exact Solutions for the Unforced Quintic Nonlinear Oscillator. Advances in Mathematical Physics, 2017, Article ID: 7396063. https://doi.org/10.1155/2017/7396063
[45]	MacColl, L.A. (1957) Theory of the Relativistic Oscillator. American Journal of Physics, 25, 535-539. https://doi.org/10.1119/1.1934543
[46]	Harvey, A.L. (1972) Relativistic Harmonic Oscillator. Physical Review D, 6, 1474-1476. https://doi.org/10.1103/PhysRevD.6.1474
[47]	Parker, E. (2017) A Relativistic Gravity Train. General Relativity and Gravitation, 49, Article No. 106. https://doi.org/10.1007/s10714-017-2267-y
[48]	Zarmi, Y. (2022) Boundary-Layer Features in Conservative Nonlinear Oscillations. Nonlinear Dynamics, 109, 2851-2863. https://doi.org/10.1007/s11071-022-07599-w
[49]	Kevorkian, J. and Cole, J.D. (1981) Perturbation Methods in Applied Mathematics. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4757-4213-8
[50]	Bender, C.M. and Orszag, S.A. (1999) Advanced Mathematical Methods for Scientists and Engineers. Springer-Verlag, Berlin. https://doi.org/10.1007/978-1-4757-3069-2
[51]	Moreau, W., Easther, R. and Neutze, R. (1994) Relativistic (An)harmonic Oscillator. American Journal of Physics, 62, 531-535. https://doi.org/10.1119/1.17513
[52]	Huang, Y.S. (1999) An Alternative Derivation of the Equations of Motion of the Relativistic Anharmonic Oscillator. American Journal of Physics, 67, 142-145. https://doi.org/10.1119/1.19209
[53]	de la Fuente, D. and Torres, P.J. (2017) A New Relativistic Extension of the Harmonic Oscillator Satisfying an Isochronicity Principle. Qualitative Theory of Dynamical Systems, 16, 579-589. https://doi.org/10.1007/s12346-016-0207-y
[54]	Gradshteyn, I.S. and Ryzhik, I.M. (2007) Tables of Integrals, Series and Products. 7th Edition, Elsevier, Amsterdam.
[55]	Fujiwara, K.M., et al. (2018) Experimental Realization of a Relativistic Harmonic Oscillator. New Journal of Physics, 20, Article ID: 063027. https://doi.org/10.1088/1367-2630/aacb5a