Exact Solutions of Nonlinear Second-Order Autonomous Ordinary Differential Equations: Application to Mechanical Systems

: Many physical processes can be described via nonlinear second-order ordinary differential equations and so, exact solutions to these equations are of interest as, aside from their accuracy, they may reveal beforehand key properties of the system’s response. This work presents a method for computing exact solutions of second-order nonlinear autonomous undamped ordinary differential equations. The solutions are divided into nine cases, each depending on the initial conditions and the system’s ﬁrst integral. The exact solutions are constructed via a suitable parametrization of the unknown function into a class of functions capable of representing its behavior. The solution is shown to exist and be well-deﬁned in all cases for a general nonlinear form of the differential equation. Practical properties of the solution, such as its period, time to reach an extreme value or long-term behavior, are obtained without the need of computing the solution in advance. Illustrative examples considering different types of nonlinearity present in classical physical systems are used to further validate the obtained exact solutions.


Introduction
The present work is concerned with the exact solutions of second-order nonlinear autonomous undamped ordinary differential equations. Many systems in nature can be modeled using these types of differential equations and so, their solutions are of immediate interest. Examples of applications include the vibrations of nonlinear mechanical systems [1], the oscillations of current or voltage in a nonlinear electric circuit [2] and the orbits of a two-body gravitational systems [3], among others.
Many researchers have investigated exact solutions for nonlinear ordinary and partial differential equations. Zhang [4] constructed exact polynomial solutions of linear second-order differential equations. Ghaemi et al. [5] studied the Hyers-Ulam stability of exact second-order linear differential equations. Lenci [1] developed exact solutions for coupled Duffing oscillators. Liao and Tan [6] discussed a general approach to obtain series solutions of nonlinear differential equations. Bartha et al. [7] computed the stable periodic orbits for the Mackey-Glass equation. Burmasheva and Prosviryakov [8] developed exact solutions to Navier-Stokes equations describing a gradient non-uniform unidirectional vertical vortex fluid flow. Carravetta [9] developed a power series solution for nonlinear ordinary differential equations having a quadratic form. Liu and Geng [10] obtained exact solutions to the systems of carbon nanotubes conveying fluid via symmetry reductions. Korman and Li [11] studied the multiplicity of positive solutions for differential equations with concave-convex and convex-concave nonlinear behavior. Cheviakov [12] developed exact closed-form solutions of a fully nonlinear asymptotic two-fluid model. Tur et al. [13] studied vortex structures with complex points singularities in two-dimensional Euler equations. Volpert et al. [14] developed exact solutions in front propagation problems with superdiffusion. MacNeil et al. [15] developed exact and approximate solutions for optical solitary waves in nematic liquid crystals. Zubarev and Zubareva [16] developed exact solutions for the shape of a two-dimensional conducting liquid drop in a non-uniform electric field. Santana et al. [17] developed closed-form solutions for the symmetric nonlinear free oscillations of pyramidal trusses.
The main idea behind the method developed in the present work consists of identifying the key properties of the solution and, from them, building a suitable parametrization that encapsulates the solution's behavior. The parametrization of the unknown function into an auxiliary function is then proven to be valid. The invertible function connecting the auxiliary function and the independent parameter is obtained, leading to the exact solution of the problem. Multiple reliable numerical methods for computing the solution of the type of differential equations considered in the present work exist in the literature [18,19]. Also, exact solutions for particular second-order ordinary differential equations are available [4,17]. However, when dealing with the general exact solution of such problems, there is a gap in knowledge that the present contribution aims to fill. The main originalities of the present contribution are listed below: • Exact solution of general second-order nonlinear autonomous undamped differential equations; • Identification of the system class from the initial conditions; • Identification of the system's general properties (e.g., period, time to reach extremes, long term behavior) before computing the solution.
The manuscript is organized as follows. The type of ordinary differential equation considered in the present work is discussed in Section 2. The solutions are divided into nine cases, depending on the ability of the system to absorb the energy provided by the initial conditions. The exact solutions are developed in Sections 2.2.1-2.2.9. The algorithm allowing for the implementation of the proposed exact solution is discussed in Section 3. Multiple examples of classical nonlinear mechanical systems are used in Section 4 to validate the developed exact solution. Finally, some conclusions and suggestions for further developments are presented in Section 5.

Exact Solutions
In this section, the initial value problem (IVP) of interest is defined, and some useful properties of the solution are discussed. The goal of the present work is to compute a function x ∈ C 2 : R → I ⊂ R satisfying the following second-order autonomous nonlinear ordinary differential equation (ODE) for a given function f ∈ C 0 : R → R: Throughout the paper, C k represents the set of functions that are k times continuously differentiable. The function x must also satisfy the following initial conditions for given {t 0 , x 0 ,ẋ 0 } ⊂ R: The second-order nonlinear ODE in Equation (1) can be classified as undamped since the function f (x(t)) does not depend directly on the rate of changeẋ(t). It can also be classified as autonomous since the function f (x(t)) does not depend directly on the independent parameter t, but only via the dependent unknown function x(t). A function x that satisfies Equation (1) together with the initial conditions in Equations (2) and (3) is said to be a solution of the initial value problem (IVP). The uniqueness of the solution is not investigated in the present contribution but can be studied via the classical Cauchy-Lipschitz theorem [20]. Also, under the condition of f ∈ C 0 , a solution is shown to always exist. Let U : R → R be defined as follows: When the differential equation in Equation (1) is used to model a physical system, the functions f and U represent the internal force and the internal energy of the system, respectively. Multiplying Equation (1) byẋ(t) and integrating over the interval (t 0 , t) providesẋ The principle of the method developed in the present work is to set the solution x as a composition of two functions z : R → I and h : R → R, that is, x(t) = z(h(t)). As it will be shown, the function z captures the qualitative behavior of the solution (Section 2.2), while the function h maps the parameter t from the domain of x to the domain of z (Section 2.1). The following conditions are imposed over the function z:

Parameter Mapping
For a fixed function z and a fixed initial parameter p 0 , the initial condition x(t 0 ) = z(h(t 0 )) = z(p 0 ) = x 0 allows setting h(t 0 ) = p 0 . Also, from the condition in Equation (8), the function g : R → R + in Equation (9) is well-defined: Replacing the composition x(t) = z(h(t)) in Equation (5) and imposing that h must be strictly increasing, thus we have The condition of h being strictly increasing could be replaced by h being strictly decreasing without any loss of generality. It is important, however, that h is strictly monotonic, so that a one-to-one relation between the domains of x and z can be later established. The equation above represents a first-order nonlinear ODE on h with initial condition h(t 0 ) = p 0 . Since g is strictly positive, an implicit solution to this equation can be obtained as follows: The relation in Equation (11) provides a unique value of h(t), ∀t ∈ R and, therefore, the function h(t) is well-defined and surjective. The implicit character of Equation (11) is reminiscent of the nonlinearity of Equation (1) and is similar to the one appearing in the Jacobi amplitude function, used to solve the nonlinear pendulum ODE [21]. Theorem 1. Fixed a function z : R → I ⊂ R satisfying the conditions in Equations (6)-(8) and defining h : R → R via Equation (11), the function x : R → R with x(t) = z(h(t)) is a solution of the IVP given in Equations (1)-(3). (2) can be verified as

Proof. The initial condition in Equation
The first derivative of x is given bẏ The initial condition in Equation (3) can then be verified aṡ Finally, the ODE in Equation (1) can be verified taking the second derivative of x:

Qualitative Behavior of the Solution
In this section, classes {z} are identified providing functions z that represent the qualitative behavior of the solution and satisfy the conditions in Equations (6)- (8). Let the limit sets S a and S b be defined as follows: From the definitions in Equations (22) and (23), U(x) =ẋ 2 0 /2 if, and only if, x ∈ S a ∪ S b . Also, S a ∩ S b = {x 0 } if, and only if,ẋ 0 = 0 and f (x 0 ) = 0, otherwise S a ∩ S b = ∅. When S a = ∅, we set x a = max(x 0 ) and, when S b = ∅, we set x b = min(S b ). The following theorem shows how the signs of f (x a ) and f (x b ) are restrained.
From the definitions in Equations (22) and (23) and the result in Theorem 2, the IVP is related to one of the nine cases presented in Table 1. Table 1. Nine possible cases of the IVP.

Case
Lower Bound Upper Bound ∃x a and f (x a ) = 0 ∃x a and f (x a ) = 0 ∃x b and f (x b ) = 0 As it will be shown, each of these cases have a particular behavior that motivates the construction of the developed exact solution. More precisely, a class of functions {z n } having a certain set of properties is defined for the nth case. An arbitrary function z n in this class does not represent an exact solution to the IVP, but encapsulates the general behavior of the solution. The link between the function z n and the actual solution of the IVP is then established via the mapping function h n , representing the evolution of the parametrization.
Due to the criteria of f (x a ) = 0 and/or f (x b ) = 0 in cases 5 to 9, the solution is sensible to the initial conditions in these cases, given that any variation on them would push the solution method to one of the other cases. Therefore, exact solutions are of particular interest since the numerical errors present in approximate methods can easily misrepresent the behavior of the solution, as will be shown in Section 4.

Case 1
Case 1 is described by S a = S b = ∅, Notice thatẋ 0 = 0 in this case, sinceẋ 0 = 0 would imply x 0 ∈ S a ∪ S b . When the IVP is used to model a physical system, the internal energy is never able to fully absorb the energy level and, hence, the system continues to evolve without bounds. These arguments motivate the following conditions over the function lim p→+∞ẋ 0 z(p) = +∞ (26) As an example of such a function, one has For the example function z in Equation (27), a unique p 0 satisfying Equation (7) can be obtained as Under the conditions of case 1, a function z ∈ C 1 : R → R that satisfies Equations (24)-(26) is bijective and also satisfies Equations (6)- (8).

Case 2
When the IVP is used to model a physical system, its internal energy is able to fully absorb the energy level only at one point to the right of the initial condition and then bounces to negative infinity. These arguments motivate the following conditions over the function z ∈ C 2 : R → (−∞, x b ]: As an example of such a function, one has For the example function z in Equation (31), a unique p 0 satisfying Equation (7) can be obtained as follows: Proof. The co-domain and continuity of z imply that z (p b ) = 0, since z(p b ) = x b is a local (and global) maximum. Therefore, Equation (29) implies that z is strictly increasing in (−∞, p b ] and strictly decreasing in [p b , +∞). Hence, z(p) = x b or z (p) = 0 imply p = p b . From the definition of case 2, U(z(p)) =ẋ 2 0 /2 implies z(p) = x b . The condition in Equation (6) is then satisfied.
Since z (p) does not approach 0 as p tends to ±∞, z is not lower bounded neither in (−∞, p b ] nor in [p b , +∞). Therefore, the restrictions z 1 : The following p 0 satisfies then the conditions in Equation (7): Applying L'Hôpital's rule, the limit as p → p b can be determined as The condition in Equation (8) is then also satisfied.
The time t b to reach the maximum value can then be determined from Equation (11) as follows:

Case 3
Case 3 is described by S a = ∅, S b = ∅ and f (x a ) < 0. When the IVP is used to model a physical system, its internal energy is able to fully absorb the energy level only at one point to the left of the initial condition and then bounces to positive infinity. It is worth noting that this case corresponds to the symmetric inverse of case 2 in Section 2.2.2 and, hence, the solutions have similar behavior. These arguments motivate the following conditions over the function z ∈ C 2 : R → [x a , +∞): As an example of such a function, one has For the example function z in Equation (38), a unique p 0 satisfying Equation (7) can be obtained as follows: Theorem 5. Under the conditions of case 3, a function z ∈ C 2 : R → [x a , +∞) that satisfies Equations (36) and (37) is surjective and also satisfies Equations (6)- (8).
Proof. The co-domain and continuity of z imply that z (p a ) = 0, since z(p a ) = x a is a local (and global) minimum. Therefore, Equation (36) implies that z is strictly decreasing in (−∞, p a ] and strictly increasing in [p a , +∞). Hence, z(p) = x a or z (p) = 0 imply p = p a . From the definition of case 3, U(z(p)) =ẋ 2 0 /2 implies z(p) = x a . The condition in Equation (6) is then satisfied.
Since z (p) does not approach 0 as p tends to ±∞, z is not upper bounded neither in (−∞, p a ] nor in [p a , +∞). Therefore, the restrictions z 1 : (−∞, p a ] → [x a , +∞) and z 2 : [p a , +∞) → [x a , +∞) are bijective. The following p 0 satisfies then the conditions in Equation (7): Applying L'Hôpital's rule, the limit as p → p a can be determined as follows: The condition in Equation (8) is then also satisfied.
The time t a to reach the minimum value can then be determined from Equation (11) as follows: The solution's behavior is lower and upper bounded. When the IVP is used to model a physical system, its internal energy is able to fully absorb the energy level at points to the left and right of the initial condition and, hence, bounces between them. These arguments motivate the following conditions over the function In Equations (44)-(46), p is the period of z. As an example of such a function, one has For the example function z in Equation (31), a p 0 satisfying Equation (7) can be obtained as follows: (44) and (45) then imply that p = p a + np or p = p b + np, with n ∈ N. The condition in Equation (6) is then satisfied.
It is assumed, without loss of generality, that p a < p b . From Equation (46) and recalling that z ∈ C 2 , z (p) > 0 in (p a , p b ) and z (p) < 0 in (p b , p a + p), the restrictions ] are then bijective. The following p 0 satisfies then the conditions in Equation (7): The limit as p → p a + np, with n ∈ N, is determined via L'Hôpital's rule, taking into account that z (p a ) > 0 (Equation (47)): The limit as p → p b + np, with n ∈ N, is determined via L'Hôpital's rule, taking into account that z (p b ) < 0 (Equation (47)): The condition in Equation (8) is then also satisfied.
The periodicity of z implies the periodicity of z and g. From Equation (11), the following relation can then be obtained: The periodicity of the solution x, with period t = h −1 (p) − t 0 , can then be shown: When the IVP is used to model a physical system, its internal energy is never able to fully absorb the energy level, but tends to/from a point (located to the right of the initial condition) and bounces from/to infinity. These arguments motivate the following conditions over the function z ∈ C 2 : R → (−∞, x b ): As an example of such a function, one has For the example function z in Equation (57), a unique p 0 satisfying Equation (7) can be obtained as follows: Under the conditions of case 5, a function z ∈ C 2 : R → (−∞, x b ) that satisfies Equations (55) and (56) is bijective and also satisfies Equations (6)-(8).
Proof. The limit in Equation (56) implies that z (ẋ 0 p) → 0 as p → +∞. Therefore, from Equation (55),ẋ 0 z (p) > 0, ∀p ∈ R and z is injective. Recalling the definition of case 5 and the co-domain of z, z(p) / ∈ S a ∪ S b , ∀p ∈ R and the condition in Equation (6) is identically satisfied. Equation (55) implies that z (ẋ 0 p) does not approach 0 as p → −∞ and hence, z is not lower bounded. The continuity of z then implies its surjectivity and so the condition in Equation (7) is satisfied. Finally, since U ∈ C 0 , z ∈ C 2 and z (p) = 0, ∀p ∈ R, the condition in Equation (8) is also satisfied.
2.2.6. Case 6 Case 6 is described by S a = ∅, S b = ∅ and f (x a ) = 0. Notice thatẋ 0 = 0 in this case, sinceẋ 0 = 0 would imply x 0 ∈ S b . When the IVP is used to model a physical system, its internal energy is never able to fully absorb the energy level, but tends to/from a point (located to the left of the initial condition) and bounces from/to infinity. It is worth noting that this case corresponds to the symmetric inverse of case 5 in Section 2.2.5 and, hence, the solutions have a similar behavior. These arguments motivate the following conditions over the function z ∈ C 2 : R → (x a , +∞): As an example of such a function, one has z(p) = x a + exp(+ẋ 0 p) For the example function z in Equation (61), a unique p 0 satisfying Equation (7) can be obtained as follows: Under the conditions of case 6, a function z ∈ C 2 : R → (x a , +∞) that satisfies Equations (59) and (60) is bijective and also satisfies Equations (6)-(8).
Proof. The limit in Equation (60) implies that z (ẋ 0 p) → 0 as p → +∞. Therefore, from Equation (55),ẋ 0 z (p) > 0, ∀p ∈ R and z is injective. Recalling the definition of case 5 and the co-domain of z, z(p) / ∈ S a ∪ S b , ∀p ∈ R and the condition in Equation (6) is identically satisfied. Equation (55) implies that z (ẋ 0 p) does not approach 0 as p → −∞ and hence, z is not lower bounded. The continuity of z then implies its surjectivity and so the condition in Equation (7) is satisfied. Finally, since U ∈ C 0 , z ∈ C 2 and z (p) = 0, ∀p ∈ R, the condition in Equation (8) is also satisfied.

Case 7
Case 7 is described by S a = ∅, S b = ∅, f (x a ) < 0 and f (x b ) = 0. When the IVP is used to model a physical system, its internal energy is able to fully absorb the energy level at a point to the left of the initial condition and tends from and to a point to the right of the initial condition that is also able to fully absorb the energy level. These arguments motivate the following conditions over the function As an example of such a function, one has For the example function z in Equation (67), a unique p 0 satisfying Equation (7) can be obtained as follows: Proof. From Equations (63)-(65) and the co-domain of z, z (p) = 0 implies p = p a and z (p) < 0 in (−∞, p a ) and z (p) > 0 in (p a , +∞). Recalling the definition of case 7, U(z(p)) =ẋ 2 0 /2 implies z(p) = x a , implying p = p a . The condition in Equation (6) is then satisfied.
The limit conditions in Equations (63) and (64) and the continuity of z imply its surjectivity in (−∞, p a ) and in (p a , +∞). Hence, the restrictions z 1 : (−∞, p a ] → [x a , x b ) and z 2 : [p a , +∞) → [x a , x b ) are bijective. The following p 0 satisfies then the conditions in Equation (7): The limit as p → p a is determined via L'Hôpital's rule: The condition in Equation (8) is then also satisfied.
From Equation (65), x(t a ) = z(h(t a )) = x a implies h(t a ) = p a . The time t a to reach the minimum value can then be determined from Equation (11) as follows: 2.2.8. Case 8 When the IVP is used to model a physical system, its internal energy is able to fully absorb the energy level at a point to the right of the initial condition and tends from and to a point to the left of the initial condition that is also able to fully absorb the energy level. It is worth noting that this case corresponds to the symmetric inverse of case 7 in Section 2.2.7 and, hence, the solution have a similar behavior. These arguments motivate the following conditions over the function z ∈ C 2 : R → (x a , x b ]: As an example of such a function, one has For the example function z in Equation (76), a unique p 0 satisfying Equation (7) can be obtained as follows:  (7): The limit as p → p b is determined via L'Hôpital's rule: The condition in Equation (8) is then also satisfied.
The time t b to reach the maximum value can then be determined from Equation (11) as follows: 2.2.9. Case 9 Case 9 is described by S a = ∅, S b = ∅, f (x a ) = 0 and f (x b ) = 0. When the IVP is used to model a physical system, its internal energy is never able to fully absorb the energy level, but tends from and to points to the left and right to the initial condition that are able to fully absorb the energy level. These arguments motivate the following conditions over the function z ∈ C 1 : R → (x a , x b ): x 0 z (p) > 0, ∀p ∈ R (83) As an example of such a function, one has For the example function z in Equation (84), a unique p 0 satisfying Equation (7) can be obtained as follows: Theorem 11. Under the conditions of case 9, a function z ∈ C 1 : R → (x a , x b ) that satisfies Equations (81)-(83) also satisfies Equations (6)-(8).

Algorithm
In this section, an algorithm to compute the exact solution of the IVP is described. The algorithm takes as inputs the function f ∈ C 0 , the initial condition parameters t 0 , x 0 andẋ 0 and the current parameter t and outputs the exact solution x(t).
The first step consists of determining the IVP type (Table 1), which is determined by checking the existence of points to left (x a ) and to the right (x b ) of x 0 such that U(x) =ẋ 2 0 /2. Unfortunately, there are not analytical methods to determine the existence of roots of a general function and so, a numerical approach is used. More precisely, the algorithm takes as additional input the boundaries (x a < x 0 and x b > x 0 ) and divides the intervals [x a , x 0 ] and [x 0 , x b ] in a uniform mesh. Next, the bisection method is used to determine if a point x a of S a exists in each sub-interval of [x a , x 0 ] and if a point Taking the direction starting with x 0 in each case ensures that, if the points x a and/or x b are found, they represent the maximum and minimum of S a and S b , respectively. The bisection method is also used in each of these sub-intervals to search for points x with f (x) = 0. If this later condition is satisfied, the condition U(x) =ẋ 2 0 /2 is checked and the special points of cases 5 to 9 can be computed.
The use of the bisection method assumes that the search function is monotonic in each of the sub-intervals. This hypothesis can be assured by controlling the size of the uniform mesh. It is worthy noting that, for specific functions U or f analytical tools can often be applied to determine the existence of the roots, excluding the need for numerical root-finding methods.
Once the type of the IVP is determined, a function z respecting the conditions in Equations (6)-(8) can be chosen. In the present work, the example functions z and initial parameters p 0 in Sections 2.2.1-2.2.9 are used; however, any function in the given class could be adopted. For a fixed z, the value of h(t) can be determined from Equation (11) and the exact solution can be computed as x(t) = z(h(t)). Due to the careful definition of the function g in Equation (9), the solution can be computed without any singularities or numerical instabilities. Moreover, since the function h is strictly monotonic, the implicit relation in Equation (11) can always be solved to an arbitrary level of precision. It is worth noting that, the use of a numerical approach is used to determine the type of the IVP and compute the mapping h(t) for a given parameter t, does not spoil the exactness of the solution x(t) since the latter can be computed to an arbitrary level of precision and no approximation of the ODE is employed, yielding a solution that is independent of the time step size.

Examples and Validation
In this section, IVPs used to model well-known nonlinear mechanical systems are used to showcase the accuracy and capacities of the developed method. The explicit fourth-order Runge-Kutta [22] and implicit Newmark [18] numerical methods are used for comparison, where a converged time step is adopted. Special initial conditions, close or equal to critical points, are used to highlight the exactness of the developed solution even when numerical methods fail due to drift errors.

Nonlinear Pendulum
The first example studied in this section is the nonlinear pendulum ( Figure 1). The pendulum consists of a point mass m attached to a rigid string with length L. The point mass is subjected to the tension force of the string and the action of gravity with acceleration g. The equation of motion of the system can be written as follows [23,24]: with initial conditions The force f and energy U functions of the system can then be obtained as follows: The internal energy function U of the nonlinear pendulum is shown in Figure 2. Three initial conditions, corresponding to cases 1, 4 and 9 from Sections 2.2.1, 2.2.4 and 2.2.9, respectively, are considered in order to highlight the capabilities of the developed exact solution in extreme cases. Considering the downward initial position θ 0 = 0, the critical initial velocityθ * 0 that brings the pendulum to the upward position (θ(t * ) = π) can be obtained from Equation (90) The first initial condition studied is taken as θ 0 = 0 andθ 0 = 1.005θ * 0 , corresponding to an energy level just above the critical one and is represented by the red line in Figure 2. The system is then never capable of fully absorbing the energy level and a certain amount of kinetic energy is always present, corresponding to the behavior described in Section 2.2.1. The results obtained with the developed exact solution, the Newmark and RK4 methods, are shown in Figure 3. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general monotonic behavior. The second initial condition studied is taken as θ 0 = 0 andθ 0 = 0.995θ * 0 , corresponding to an energy level just below the critical one and is represented by the blue line in Figure 2. The system is then capable of fully absorbing the energy level at points to the left and to the right of the initial position, corresponding to the behavior described in the Section 2.2.4. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 4. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general oscillatory behavior. It is worth noting that, while both numerical methods starts to diverge from the periodic behavior, the exact solution retains this key property. Moreover, the dashed black lines in Figure 4 represent the period t = 4.748945 s, calculated via Equation (54), showcasing its accuracy.  Although the two numerical methods appear to drift by the same amount, in Figure 4 the red (RK4) and black (Newmark) points vary for the same values of the parameter t, and for t > 6 these points start to visibly differ from one another, indicating loss of accuracy of the numerical methods. Moreover, the accuracy of the exact solution can be verified by the dash lines indicating the period of the response. Since the energy provided by the initial conditions in this case is smaller than the critical one, the solution must always be periodic. Until the second dashed line (t ≈ 9.5), all solutions have approximately the same period. However, after this point, the solutions from the two numerical methods (via drift errors) have a smaller period, while the response from the exact solution maintain the original period.
For a null initial velocity (θ 0 = 0), the response of the system is always periodic, falling in case 4 of Section 2.2.4. The period of the response can than be calculated as a function of initial positions θ 0 as shown in Figure 5. These calculations are performed without the need of computing the time response of the system in advance and highlight an advantage of exact methods. The third initial condition studied is taken as θ 0 = 0 andθ 0 =θ * 0 , corresponding to an energy level equal to the critical one and is represented by the green line in Figure 2. The system is then never capable of fully absorbing the energy level but approaches such points to the left and to the right of the initial position, corresponding to the behavior described in the Section 2.2.9. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 6. As can be observed, the developed exact solution shows a good agreement with the numerical methods at the beginning. However, as this initial condition corresponds to a critical one, even small numerical drift errors will throw the numerical response to either case 1 or case 4 described above. These phenomena can be clearly observed in Figure 6 as the Newmark method drift the solution to case 1 (monotonic behavior), while the RK4 method drifts the solution to case 4 (periodic behavior). The developed exact solution, however, retains the asymptotic behavior of the system, highlighting once more the advantages of exact methods.

Nonlinear Pyramidal Truss
The next example studied in this section is the nonlinear pyramidal truss (Figure 7). The truss consists of n bars attached to central point, at a height H, that can only move vertically (z) and are blocked in the base, that forms a circle with radius R. The central point is subjected to the tension force of all the bars. The bars have a section with area A and are made of a material with elastic modulus E and specific mass ρ. The natural frequency ω 0 of the system in the reference configuration is then given by [17] where L = √ R 2 + H 2 . Considering a logarithmic strain measure and the adimensional time t = ω 0 t and position z = z/H, the equation of motion of the system can be written as follows [17]: where α = R/H and β = L/H = (a) (b) The force f and energy U functions of the system can then be obtained: The internal energy function U of the nonlinear pyramidal truss is shown in Figure 8. Two initial conditions, corresponding to cases 7 and 8 from Sections 2.2.7 and 2.2.8, respectively, are considered in order to highlight the capabilities of the developed exact solution in extreme cases. The truss has three static equilibrium points, that is, zeros of the force function f . The equilibrium points z = ±1 are stable, while the equilibrium point z = 0 is unstable, as indicated by the curvature of this points in Figure 8. Starting from any of the stable equilibrium points (z = ±1) the critical velocityż * 0 required to reach the unstable equilibrium point (z = 0) can be obtained: The first initial condition studied is taken as z 0 = −1 andż 0 = +ż * 0 , corresponding to an energy level equal to the critical one and is represented by the blue line in Figure 8. The system is capable of fully absorbing the energy level in a point to left of the initial condition and approaches such a point to the right of the initial condition, corresponding to the behavior described in the Section 2.2.7. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 9. As can be observed, the developed exact solution shows a good agreement with the numerical methods at the beginning. However, as this initial condition corresponds to a critical one, even small numerical drift errors will throw the numerical response to case 4, either oscillating around just one stable equilibrium point or around both of them. These phenomena can be clearly observed in Figure 9 as the Newmark method drifts the solution to case 4. For a larger time interval, the RK4 also showed such drift numerical errors. The developed exact solution, however, retains the asymptotic behavior of the system, highlighting once more the advantages of exact methods. Moreover, the dashed black line in Figure 9 represents the time t = 1.696756 s to reach the minimum value x a = −1.732043, calculated via Equation (71), showcasing its accuracy. The second initial condition studied is taken as z 0 = +1 andż 0 = +ż * 0 , corresponding to an energy level equal to the critical one and is represented by the red line in Figure 8. The system is capable of fully absorbing the energy level in a point to right of the initial condition and approaches such a point to the left of the initial condition, corresponding to the behavior described in Section 2.2.8. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 10. As can be observed, the developed exact solution shows a good agreement with the numerical methods at the beginning. However, as this initial condition corresponds to a critical one, even small numerical drift errors will throw the numerical response to case 4, either oscillating around just one stable equilibrium point or around both of them. These phenomena can be clearly observed in Figure 10 as the Newmark method drifts the solution to case 4. For a larger time interval, the RK4 also showed such drift numerical errors. The developed exact solution, however, retains the asymptotic behavior of the system, highlighting once more the advantages of exact methods. Moreover, the dashed black line in Figure 10 represents the time t = 1.696756 s to reach the maximum value x b = 1.732043, calculated via Equation (80), showcasing its accuracy.

Spike System
The next example of this section consists of a system where the internal energy spikes around an unstable equilibrium point. This kind of behavior is common in chemical systems or probabilistic models [25][26][27]. The equation of motion of the system can be described as follows:ẍ The force f and energy U functions of the system can then be obtained: The internal energy function U of the spike system is shown in Figure 11. Four initial conditions, corresponding to cases 2, 3, 5 and 6 from Sections 2.2.2, 2.2.3, 2.2.5 and 2.2.6, respectively, are considered in order to highlight the capabilities of the developed exact solution in extreme cases. The system has one unstable equilibrium point at x = 0.
The first initial condition studied is taken as x 0 = −2 andẋ 0 = +2 2U(0.5), corresponding to an energy level that makes the system reaches the maximum point x b = −0.50 and is represented by the blue line in Figure 11. The system is capable of fully absorbing the energy level only to a point to right of the initial condition, corresponding to the behavior described in the Section 2.2.2. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 12. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general parabolic behavior. Moreover, the dashed black line in Figure 12 represents the time t = 0.2707428 s to reach the maximum value, calculated via Equation (35), showcasing its accuracy.
The second initial condition studied is taken as x 0 = +2 andẋ 0 = −2 2U(0.5), corresponding to an energy level that makes the system reach the minimum point x b = +0.50 and is represented by the red line in Figure 11. The system is capable of fully absorbing the energy level only to a point to left of the initial condition, corresponding to the behavior described in Section 2.2.3. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 13. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general parabolic behavior. Moreover, the dashed black line in Figure 13 represents the time t = 0.2707428 s to reach the maximum value, calculated via Equation (42), showcasing its accuracy.
The third initial condition studied is taken as x 0 = −2 andẋ 0 = +2 2U(0), corresponding to an energy level equal to the critical one, making the system approach the critical point x = 0 and is represented by the green line in Figure 11. The system is never capable of fully absorbing the energy level but approaches such a point to right of the initial condition, corresponding to the behavior described in the Section 2.2.5. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 14. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general asymptotic behavior.
The fourth initial condition studied is taken as x 0 = +2 andẋ 0 = −2 2U(0), corresponding to an energy level equal to the critical one, making the system approach the critical point x = 0 and is represented by the magenta line in Figure 11. The system is never capable of fully absorbing the energy level but approaches such a point to left of the initial condition, corresponding to the behavior described in Section 2.2.6. The results obtained with the developed exact solution, the Newmark and RK4 methods are shown in Figure 15. As can be observed, the developed exact solution shows a good agreement with the numerical methods and represents well the general asymptotic behavior.

Discussion and Conclusions
In this work, exact solutions regarding second-order nonlinear autonomous undamped differential equations were developed. The behavior of the solution was linked to its energy level and nine classes of solutions were identified. The exact solutions were constructed via a suitable parametrization of the unknown function to a capable class that encapsulates its nonlinear response. The solution was shown to exist and be well-defined over all cases, for a general nonlinear form of the differential equation. Illustrative examples considering different types of nonlinearity were used to further validate the obtained exact solutions.
In future works, the exact solution could perhaps be extended to damped systems. For undamped systems, the energy level is constant and so, the topology of the solution remains the same as the system evolves. For damped systems, this is not generally the case. The solution can pass from one topological form to another as the energy level of the system changes. These considerations could be taken into account when proposing a new parametrization to the solution.
The extension of the exact solution to non-autonomous systems or systems with multiple degrees of freedom is not trivial. In these cases, the solution can be chaotic and so, identifying a general behavior to propose a parametrization seems to be a very difficult task. For the non-autonomous case, a scalable form could perhaps be used to parametrize the solution. In this case, attention should be made to the change of suitable form of the parametrization as the solution evolves.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: IVP initial value problem ODE ordinary differential equation