Predicting saddle-node bifurcations using transient dynamics: a model-free approach

This paper proposes a novel method for predicting the presence of saddle-node bifurcations in dynamical systems. The method exploits the effect that saddle-node bifurcations have on transient dynamics in the surrounding phase space and parameter space, and does not require any information about the steady-state solutions associated with the bifurcation. Specifically, trajectories of a system obtained for parameters close to the saddle-node bifurcation present local minima of the logarithmic decrement trend in the vicinity of the bifurcation. By tracking the logarithmic decrement for these trajectories, the saddle-node bifurcation can be accurately predicted. The method does not strictly require any mathematical model of the system, but only a few time series, making it directly implementable for gray- and black-box models and experimental apparatus. The proposed algorithm is tested on various systems of different natures, including a single-degree-of-freedom system with nonlinear damping, the mass-on-moving-belt, a time-delayed inverted pendulum, and a pitch-and-plunge wing profile. Benefits, limitations, and future perspectives of the method are also discussed. The proposed method has potential applications in various fields, such as engineering, physics, and biology, where the identification of saddle-node bifurcations is crucial for understanding and controlling complex systems.


Introduction
Stability is a critical property in engineering, indicating a system's ability to return to its original state after experiencing an arbitrarily small disturbance.Typically, stability is used to guarantee the safety of a system in a dynamical sense, but this is not a sufficient condition.When a system operates outside the region in which it can be considered linear, trajectories may not converge back to the stable solution.Instead, they may lead to another solution, which can be detrimental to the system, or even diverge completely.This could potentially compromise the physical integrity of the system and lead to an accident.
Subcritical Andronov-Hopf bifurcations provide a typical scenario leading to equilibrium points only locally stable.Let us consider the bifurcation diagram depicted in Fig. 1.The branch of unstable periodic solutions emerging at the bifurcation coexists with the locally stable trivial solution.Then, it turns to the right in correspondence of a saddle-node (or fold) bifurcation.After the saddle-node bifurcation, the branch becomes stable and reaches the region where the trivial solution is unstable.The bifurcation parameter space of a system exhibiting a similar bifurcation diagram can be divided in three regions, marked as A, B and C in the figure.Let us assume that the trivial solution is desirable, while any other steady-state solution is dangerous for the system integrity.In region A, the trivial solution is globally stable (assuming that no other solution exists).In this region, if the system is in equilibrium and it is subject to an external perturbation, it will return to the equilibrium state, regardless of the intensity of the perturbation.From a dynamical perspective, we can say that this region is safe.In region C, the trivial solution is unstable.Therefore, the system will never remain in equilibrium there, but it will converge towards the stable periodic solution, whatever the Fig. 1 Typical scenario of a subcritical Andronov-Hopf and saddle-node bifurcations.Solid line: stable solutions, dashed lines: unstable solutions initial conditions are.This region is obviously unsuitable for safe operation.Region B is the most dynamically involved.Two stable solutions coexist-the trivial and a periodic one-and the unstable periodic solution lies on the boundary of the basins of attraction of the two stable solutions.If the system is in the trivial solution state, a sufficiently large perturbation might make the system leave the basin of attraction of the trivial solution and converge towards the periodic solution.However, a small perturbation will only lead to transient effects, and the system will converge again to the trivial solution.Accordingly, although the system can correctly work in this region, it is unsafe as dangerous stationary periodic motions can arise.The trivial solution is only locally stable and not globally.We remark that a stable periodic solution coexisting with the trivial one might never occur for long operational time, and appear suddenly if the system reaches a specific region of the phase space, causing an accident hardly explainable.
Although typical industrial approaches may overlook region B, it can be found through a bifurcation analysis and continuation of the branch of periodic solutions, either numerically or analytically.Experimentally, or through direct numerical simulations, region B can be identified by performing a sweep-up and a sweep-down of the bifurcation parameter.However, this procedure is not always possible, as persistent periodic oscillations might undermine the integrity of the system, preventing this approach, as for instance for flutter instabilities [25].
Fig. 2 Typical isolated branch of limit cycle oscillations.Solid line: stable solutions, dashed lines: unstable solutions Figure 2 illustrates a similar, but much more subtle bifurcation scenario.In this case, the branches of periodic solutions are closed in a loop, forming a socalled isola, or detached branch [26].These are encountered in numerous systems, such as traffic models [27], shimmy of rigid wheels [28], towers and bridges under wind excitation [29,30], chemical reactors [31,32], and predator-prey ecosystems [33,34].In a similar bifurcation scenario, the trivial solution is stable for the whole considered range of the bifurcation parameter.Region A and C are both globally stable and topologically equivalent.Conversely, region B has similar properties as region B in the previous example.However, in this case, the identification of region B is much more complicated.In fact, both a sweep up or a sweep down of the bifurcation parameter will not reveal its existence, either in experiments or from a numerical model.Neither numerical continuation would provide better results.In some cases, an analytical investigation of the system's dynamics might reveal the existence of the isola [26].Otherwise, one of the few ways to identify the isola is to perform a global analysis of the system dynamics [35,36], which might be very timeconsuming.We note that usually isolas are not fully detached from other branches of system's solutions; however, their merging might be in unexplored regions of the parameter space [37].
As the two scenarios above illustrate, dangerous regions where solutions are only locally stable are often delimited by saddle-node bifurcations.However, identifying these bifurcations is not always trivial, as classical methods to find them, whether analytical, numerical [38,39], or experimental [40], require first identifying one of the branches of periodic solutions forming them.Accordingly, they are often undetected.
Several methods exist to predict the loss of stability of an equilibrium from the system's transient dynamics in the vicinity of the stability loss, which are also used in engineering practice [41].Similar techniques can also be implemented for forecasting fold bifurcations of equilibrium points [42,43]-which are usually referred to as tipping points in the context of climate [44][45][46] and ecological systems [47,48].However, they predict the fold only from the stable branch directly leading to the fold, as the classical methods discussed above do, and not from the globally stable region.
Recently, Lim and Eupreano [49] developed a method for forecasting Andronov-Hopf bifurcations and the emerging bifurcation branch only from a few trajectories obtained in the pre-bifurcation scenario.The technique leverages the so-called critical slowing down [50], i.e., the increased relaxation time of a dynamical system near a bifurcation.The method was further developed and tested on various systems, including aeroelastic flutter [51][52][53], parametrically excited systems [54] and traffic jams [55].In the case of subcritical Andronov-Hopf bifurcations, the method enables to identify also saddle-node bifurcations from the estimated bifurcation branch.Between existing methods to understand a system's dynamics from time series without resorting to system identification, it is worth mentioning the works by Xu et al. [56,57] for the localization of saddle points.
The objective of this study is to provide a method to identify saddle-node bifurcations, which does not require any knowledge about the branch of periodic solutions forming it.Instead, it exploits the effects that saddle-node bifurcations have on the surrounding phase space and parameter space, similarly to the method proposed in [49], but specifically focused on saddle-node bifurcations.As such, the method does not even require knowledge of the equations of motion of the system, but only a few trajectories in the vicinity of the bifurcation.

Fundamental idea
We consider the linear loss of stability of an equilibrium point, which occurs for a specific critical value μ cr of the parameter μ.The equilibrium is stable for μ < μ cr and unstable for μ > μ cr .The stability loss is related to the real part of an eigenvalue (or of a couple of eigenvalues) of the Jacobian matrix that becomes positive.In most cases, as μ approaches μ cr , the real part of the eigenvalues, initially negative, increases; then reaches zero for μ = μ cr , and finally becomes positive for μ > μ cr .Since the real part of this eigenvalue is strictly related to the local speed of converge to the equilibrium, the speed of convergence will decrease as μ approaches μ cr , will be nearly zero for μ = μ cr (the equilibrium is non-hyperbolic), and will be negative (trajectories diverge from the equilibrium) for μ > μ cr .If this scenario is verified, although the bifurcation separates two qualitative distinct scenarios, the variation of the speed of convergence is smooth.This feature can be used to predict stability loss in real-time by computing the speed of convergence [41].Now we consider a generic 2-dimensional system, such as a single-degree-of-freedom (DoF) mechanical system, that exhibits a saddle-node bifurcation as the one separating regions A and B in either Figs. 1 or 2. For simplicity of explanation, let us assume that, if we neglect the phase, the system can be reduced to its 1dimensional polar coordinate form given by where steady-state periodic solutions are given by equilibrium points of ρ, while the trivial equilibrium is given by ρ = ρ 0 := 0. In region A (μ < μ sn ), the only steady-state solution is the trivial solution ρ 0 that is globally stable, i.e., trajectories converge towards it from any point of the phase space (any initial ρ(0) value).Referring to Fig. 3, region A corresponds to subplots (a) and (b).
At the saddle-node bifurcation (μ = μ sn , Fig. 3c), the stable and unstable periodic orbits merge into the saddle-node (ρ sn ).Accordingly, trajectories starting with ρ(0) > ρ sn will converge towards the saddlenode point; otherwise, they will go towards the trivial solution.
In region A, oscillation amplitude monotonically decreases for any ρ.While in region B, it decreases for ρ > ρ 2 or for ρ < ρ 1 , while it increases for ρ 1 < ρ < ρ 2 .The saddle-node bifurcation marks the boundary between these two regions topologically different.Considering that, in the vicinity of the bifurcation, the system transitions from a decreasing to an increasing amplitude behavior, we hypothesize that the transition is gradual, similar to an equilibrium losing stability, as discussed above.This hypothesis implies that, for amplitudes close to the saddle-node bifurcation, the rate of the amplitude decay, measured, for instance, by the logarithmic decrement, gradually decreases as we approach the bifurcation, it is zero at the bifurcation, and it progressively becomes negative after the bifurcation, for μ > μ sn .If this hypothesis is verified, then the trend of the logarithmic decrement of trajectories converging towards a stable solution can be utilized to predict the presence of a saddle-node bifurcation.The next section develops a formulation for verifying and exploiting this phenomenon through an example.Then, the method is applied to other systems to analyze the generality and limitations of the method.

Saddle-node bifurcation prediction
We consider a single-DoF system having a nonlinear damping force, whose dynamics is governed by the following differential equation where c 1 and c 3 are real numbers.This system was thoroughly studied in [26].For c 1 positive, the trivial solution is always stable.However, for c 3 = c * 3 := 40c 1 /9 a saddle-node bifurcation exists, generating one branch of stable periodic solutions and one branch of unstable ones for c 3 > c * 3 , as illustrated in Fig. 4. In the following, we fix c 1 = 0.1, which leads to c * 3 = 0.444.For c 3 < c * 3 , the trivial solution is globally stable, while for c 3 > c * 3 , it is only locally stable.For c 3 = 0, the system corresponds to a linear damped oscillator.A free vibration decay for this case is illustrated in Fig. 5a.The system experiences a classical exponential decay.The system's trajectory in the phase space generates a logarithmic spiral (Fig. 5d), which is self-similar.We computed the logarithmic decrement as where A i is the i th amplitude peak.The logarithmic decrement is constant for a single-DoF linear system, as visible from the line marked with c 3 = 0 in Fig. 6a.We now consider the case of c 3 = 0.4, which is relatively close to the saddle-node bifurcation occurring at c 3 = c * 3 = 0.444.Figure 5b illustrates the free vibration decay and the relative trajectory in the phase space in Fig. 5e.The envelope of the displacement peaks does not produce an exponential decay anymore.In particular, the envelope presents a sort of bump for x ≈ 0.7.The spiraling trajectory in the phase space (Fig. 5e) presents denser lines for the same oscillation amplitude.The logarithmic decrement (Fig. 6a, c 3 = 0.4) is not constant anymore, but it has a clear minimum for x ≈ 0.8.This observation confirms, at least for this specific case, the hypothesis formulated in Sect.2, i.e., that the variation of the logarithmic decrement is gradual as the saddle-node bifurcation is approached and not sudden.
Considering the saddle-node bifurcation, we expect that, for c 3 ≥ c * 3 , trajectories with large initial conditions will not converge towards the trivial solution but to the stable periodic solutions.Once they reach the periodic solution, the logarithmic decrement should be zero.This scenario is confirmed by the time series, phase portraits, and logarithmic decrement curves represented in Figs.5c, f and 6a, respectively, obtained for c 3 = c * 3 .Comparing the phase portraits in Figs.5e and  f, we also note that the periodic solution lays exactly in the phase space region where the spiral in Fig. 5e is thicker.
According to the depicted scenario, the presence of a local minimum in the logarithmic decrement curve can be interpreted as a possible sign of a saddle-node bifurcation approaching.By performing several simulations of the system for increasing c 3 and plotting the locus of the minima of the logarithmic decrement with respect to the bifurcation parameter c 3 , a curve is obtained that intersects the zero level for c 3 = c * 3 , as depicted in Fig. 6b.In this specific case, the curve is practically a straight line, suggesting that the critical value of the saddle-node bifurcation can be accurately estimated by obtaining only two points of the curve.Since the logarithmic decrement curves are defined by a fairly small number of points-which increases for c 3 approaching c * 3 -their minima are obtained through a spline interpolation.
The trend of the logarithmic decrement with respect to the system energy for various c 3 values is represented by the color map in Fig. 6c.The color map clearly shows that, as the saddle-node bifurcation and the periodic solution is approached, the logarithmic decrement tends to zero.Between the two branches of periodic solutions, for c 3 > c * 3 , the logarithmic decrement is positive (not illustrated in the figure).
This example proves that the supposed scenario in the vicinity of the saddle-node bifurcation exists and can be exploited for predicting the bifurcation.However, the considered system is relatively simple, as it has only one DoF, and nonlinearities are smooth and polynomial.
In the following sections, more challenging cases will be considered-namely, a non-smooth, a timedelayed, and a multi-DoF system.

Non-smooth system: mass-on-moving-belt
We consider the classical mass-on-moving-belt system represented in Fig. 7 [11,12], which is an archetypal model for friction-induced vibrations investigation, typically used for studying brake squeal [58] and violin string dynamics [59].
The non-dimensional equations of motion of the system are [11] where ζ is the linear damping ratio and The friction coefficient μ is given by the exponential decaying function where v rel = v − ẋ is the relative velocity.The friction force presents a discontinuity for v rel = 0.For this study the same parameter values utilized in [12] are adopted, i.e., μ s = 1, μ d = 0.5, v 0 = 0.5 and ζ = 0.05.As thoroughly explained in [11], the system's equilibrium loses stability through a subcritical Andronov-Hopf bifurcation for a specific belt velocity v = v cr .
The equilibrium is stable for v > v cr and unstable otherwise.A branch of unstable periodic solution arises at v = v cr .This branch then merges with a branch of stable periodic solutions through a saddle-node bifurcation, as illustrated in Fig. 8.We denote the velocity corresponding to the saddle-node bifurcation with v * .For the adopted parameter values, v cr = 1.151 and v * = 1.83.
The system dynamics is topologically similar to the one presented in the previous section.However, in this case, the stable periodic solution is characterized by a stick-slip motion.Because of the discontinuity of the system, the branches of stable and unstable periodic solutions merge in a non-smooth way.In the following, we apply the same procedure adopted in the previous section for predicting the presence of the saddle-node bifurcation.
Figure 9 illustrates time series and trajectories in the phase space for various belt velocity v values.For v = 3.3, the saddle-node bifurcation is relatively far.From the free vibration decay (Fig. 9a) it is hard to spot any non-smoothness of the displacement curve.In the corresponding phase portrait (Fig. 9d) a small non-smoothness can be identified each time ẋ = v (v rel = 0).However, the spiraling is rather regular.
For v = 1.9, the system is much closer to the saddlenode bifurcation.The free vibration decay in Fig. 9b shows that the first three peaks are very similar to the case of v = 3.3.From Fig. 9e, we note that these three peaks occur before the system experience sticking condition.The rest of the time series in Fig. 9b is very different from the one in Fig. 9a.In fact, after the first sticking, the energy decay seems much slower, as the fourth, fifth, and sixth peaks have very similar amplitudes.The trajectory in Fig. 9e exhibits dense curves close to the first sticking point.
Finally, Figs.9c and f refer to a belt velocity v = 1.8, for which the stable and unstable periodic solutions coexist with the equilibrium.Interestingly, also in this case, the first three peaks of the free vibration decay are almost identical to the other two cases.However, as the first sticking condition is met, the system dwells in a self-sustained stick-slip periodic motion.
The phenomena observed in Fig. 9 suggest that the saddle-node bifurcation and variations of the belt velocity do not significantly affect the large amplitude dynamics (before the first sticking occurs) but mainly the small amplitude dynamics (after the first sticking).Accordingly, for predicting the saddle-node bifurca- Figure 10a shows the logarithmic decrement of the system's displacement after the first sticking for various belt velocities v.For very large velocities (v = 10), the logarithmic decrement is practically constant, as for a linear system.However, as v is reduced and the system approaches the saddle-node bifurcation, the logarithmic decrement increases as the oscillation amplitude decreases.The envelope of the minima of the logarithmic decrement curves forms a path (marked in red in the figure) that intersects the zero axis for v = v * , where the saddle-node bifurcation exists.This result implies that the same strategy that was implemented for the smooth system studied in Sect. 3 can be adopted for this discontinuous system as well.However, only the dynamics after the first sticking should be taken into account.
By collecting the minima of the logarithmic decrement curves for different belt velocity values and plotting them with respect to the belt velocity, a path is obtained, which intersects the zero level in correspondence with the saddle-node bifurcation.Accordingly, it can be used to predict for which value of v it will occur (Fig. 10b).Differently from the previous system studied in Sect.3, we note that for this system the locus of the minima is not a straight line, which makes the prediction of the exact location of the saddlenode bifurcation slightly more complicated.Nevertheless, we remark that signals that a saddle-node bifurcation might occur are already visible for velocity values very far from the critical one.For instance, comparing the logarithmic energy trend for v = 3.3 and v = 2.5 in Fig. 10a, the possible occurrence of a saddlenode bifurcation can be inferred.The amplitude of the saddle-node bifurcation is also predicted by the inter-Fig.10 The figure refers to the mass-on-moving-belt system in Eq. ( 4). a Logarithmic decrement for various v values, initiated from their minimal point (after the first sticking phase); the magenta point indicates the critical v value, the red line is the locus of the minimum of the curves.b Locus of the minima of the logarithmic decrement curves with respect to the belt velocity v section of the red curve in Fig. 10a with the zero axis.The agreement with the amplitude of the saddle-node illustrated in Fig. 8 is excellent.
Referring to the curves in Fig. 10a, we note that at low energy, the curves do not converge to the same value, as it was the case in Fig. 6a.The value of the logarithmic decrement for very small amplitudes is strictly related to the linear damping of the system.For the nonlinearly damped oscillator (Fig. 6a), this value depends only on c 1 (see Eq. 2), which was kept constant.Conversely, for the mass-on-moving-belt, the linear damping decreases as v approaches v cr ; therefore, the logarithmic decrement at small energy also decreases.For v → v cr , the logarithmic decrement at small amplitudes tends to zero, a property which can be used to predict stability loss [41].

Time-delayed system: inverted pendulum
As a third system to analyze, we consider an inverted pendulum, whose vertical position is stabilized by a torque proportional to the pendulum's position and velocity (Fig. 11).We assume that the control torque is applied with a time delay τ with respect to the measured quantities.The non-dimensional equations of motion of the system are where p and d are the proportional and differential control gains, respectively.If τ → 0, the trivial solution is asymptotically stable for p > 1 and d > 0. Conversely, for τ > 0, the trivial equilibrium becomes unstable for too large p and d values.Several studies [60,61] investigated the stability of this and similar models.In particular, it was proven that, in most cases, the system loses stability through a subcritical Andronov-Hopf bifurcation [62,63], as illustrated in Fig. 12 for d = 2 and τ = 0.5.The bifurcation diagram was obtained through a collocation method combined with pseudo-arch-length continuation via the MATLAB toolbox DDEbiftool [64].
The bifurcation scenario in Fig. 12 is qualitatively similar to the one in Fig. 1, with the difference that several saddle-node bifurcations are encountered along the branch of periodic solutions.Additionally, although this system has one DoF as the other systems studied so far, the time delay makes it infinite dimensional [65].
According to the bifurcation diagram in Fig. 12, the lowest saddle-node bifurcation occurs for p = p * = 2.126.Figure 13 depicts free vibration decay for various bifurcation parameter values p.No particular feature can be noticed in the time series for p = 1.8 (Fig. 13a).The logarithmic decrement trend shown in Fig. 14a is relatively flat, except at very low amplitudes; however, it exhibits a small indentation at ϕ ≈ 5.
Approaching the saddle-node, for p = 2.1 (Fig. 13b), the free vibration decay presents a small bump for ϕ ≈ 10 and very slow decay for ϕ ≈ 5. Figure 14a shows that, for p = 2.1, the logarithmic decrement curve presents three dips, one for each of the saddlenode bifurcations present in the bifurcation diagram for the considered oscillation amplitude.Finally, for p = 2.13 ( p > p * , Fig. 13c), the free vibration decay settles on a periodic solution; however, before settling, a small bump is visible for ϕ ≈ 10.In Fig. 14a, we notice that, for p = p * , the logarithmic decrement in correspondence with the leftmost valley goes to zero.However, also in this case, two dips are visible in correspondence with the other saddle-node bifurcations.This scenario is consistent with what was observed for the previous systems.
Collecting the leftmost local minimum of each logarithmic decrement curve and plotting them with respect to the bifurcation parameter p, we obtain an almost straight line that intersects the zero axis in correspon-dence with the saddle-node bifurcation, analogously to what was observed in the previous cases.This suggests that, also for this time-delayed system, the minimum point of the logarithmic decrement can be used to predict the occurrence of saddle-node bifurcations and its amplitude.
By computing the logarithmic decrement for a range of p values and interpolating the curves through splines, we obtained the color map shown in Fig. 14c.The minimum of the logarithmic decrement at the three depicted saddle-node bifurcations can be observed.We note a mismatch between the bifurcation diagram obtained through the DDEbiftool MATLAB toolbox [64] and the direct numerical simulations performed through the MATLAB delayed differential equation solver DDE23.The numerical simulations indicate slightly lower p values for the saddle-node bifurcations.This mismatch could be due to the relatively poor tolerance used for the integration, as tightening the tolerance significantly slows down the computations.Nevertheless, the results are qualitatively in excellent agreement.

Multi-DoF system: aeroelastic flutter
As a last example, we consider a pitch-and-plunge wing profile with an attached nonlinear tuned vibration absorber, as the one studied in [10].The pitch-andplunge model considered, used to describe an airfoil motion, was implemented in various studies [66][67][68], while the nonlinear tuned vibration absorber is practically a tuned mass damper also encompassing a nonlinear restoring force [69][70][71].Non-dimensional equations of motion governing the dynamics of the system are where y indicates the heave displacement, α the pitch rotation, and x the absorber displacement, non-dimensionalized in relation to the semichord of the airfoil, while u is the non-dimensional flow velocity.For the physical meaning of all the other parameters and a schematic mechanical model, we address the interested reader to [10]; the same nomenclature is utilized here for facilitating the comparison.The adopted parameter values are x α = 0.2, r α = 0.5, β = 0.2, ν = 0.08, = 0.5, For the considered parameter values, the bifurcation diagram for variations of the flow velocity u is shown in Fig. 15.Although the equilibrium position loses stability through a supercritical Andronov-Hopf bifurcation, the emerging branch of stable periodic solutions turns back and moves into the region where the trivial solution is stable (for u < u cr ), until it reaches a saddle-node bifurcation for u = u * .Then, the branch turns again and moves back to the unstable region.For the chosen parameter values, u cr = 1.255 and u * = 0.8205.The overall scenario in the vicinity of Fig. 15 Bifurcation diagram of the system in Eq. (8) the saddle-node bifurcation is qualitatively similar to the one represented in Figs. 1, 2, 4 or 12.However, it presents the significant difference that the system has three DoF; therefore, it is 6-dimensional.Besides, modal interactions are present between the wing and the tuned mass damper.
To analyze the effect of approaching the saddle-node bifurcation on the system's global dynamics, we perform several time integrations of the system for various u values (Fig. 16).Comparing the three time series obtained for u = 0.6 (far from the saddle-node bifurcation), u = 0.8 (very close to the bifurcation), and u = 0.822 (just after the bifurcation), we note that, for u = 0.6, the free vibration decay is quite regular, although the signal presents various vibration frequencies.For u = 0.8, the system oscillates with an α amplitude of about 0.4 for some time before converging to the equilibrium.Finally, for u = 0.822, the system settles on a stable periodic solution.The behavior is qualitatively similar to what was observed for the previously studied systems and seems promising for adopting the same procedure to estimate the saddle-node bifurcation position.
At first, we tried to reproduce the analysis performed for the previous systems with the logarithmic decrement of the system's displacement.However, as the system has more than one DoF, it is ambiguous which variable should be chosen.Besides, beating phenomena may occur, hiding the oscillation amplitude decrement if only one variable is considered.Accordingly, we decided to compute the logarithmic decrement of the system's mechanical energy, which considers the variation of the oscillation amplitude of all the system's state variables.The logarithmic energy decrement ( E ) was computed for each peak of the heave motion.
The mechanical energy of the system, excluding the contribution of aeroelastic forces, is defined as E = T + U , where the kinetic and potential energies T and U are and Figure 17a depicts the logarithmic energy decrement for various values of u.All curves present a dip that becomes increasingly pronounced as u approaches u * .Finally, the logarithmic energy decrement reaches zero as u = u * .This behavior is similar to what was observed in previous cases.However, the curves overlap and are irregular.This effect is likely related to the multi-DoF nature of the system and to the modal interaction.Although each point of the curves corresponds to a peak of the heave motion y, they are not as regular as those observed for a single DoF system.
Figure 17b shows a color map representing the logarithmic energy decrement for various energy and flow velocity u values.The depicted scenario exhibits a general trend similar to the previous examples.A clear reduction of the logarithmic energy decrement is visible near the saddle-node bifurcation and the periodic solutions, where it goes to zero.However, many vertical lines with different colors disturb the regularity of the color map.These lines are related to trajectories that do not pass near the saddle-node bifurcation.For topological reasons, this phenomenon cannot occur in a 2-dimensional system.Conversely, it is quite common in a multi-dimensional system, such as the present one.
To further analyze this effect, we consider two trajectories obtained for the same parameter value set (u = 0.79027) with slightly different initial conditions.Figure 18a and c refer to initial conditions x(0) = [2, 2, −1.9] and ẋ = 0, while Fig. 18b and d   trends.The time series in Fig. 18a shows a clear slowdown of the energy decay for α ≈ 0.4.In the phase portrait (Fig. 18c), this corresponds to denser curves in the region where the stable periodic solution will appear (marked in red in the figure).The trajectory depicted in Fig. 18d does not exhibit this feature, and the corresponding time series does not seem to show any slow-down of the energy decay.In fact, the trajectory in Fig. 18d does not transit in the vicinity of the saddle-node bifurcation.However, this is not easily visualized because of the large dimension of the system.In some sense, the trajectory in Fig. 18b and d carries weak information about the saddle-node bifur-cation, which cannot be easily captured with the technique used so far.We remark that also the technique developed in [49] would probably fail for a similar trajectory.
Regarding the prediction of the saddle-node bifurcation, Fig. 17c shows the collection of minimal points from all logarithmic energy decrement curves used for Fig. 17b (points with negative logarithmic energy decrement were excluded).Most points tend to align on a curve indicating the location of the saddle-node bifurcation.Although many points deviate significantly from this curve, introducing noise into the data and Fig. 18 Time series and phase portraits of the system in Eq. ( 8) for u = 0.79027.Initial conditions are ẋ(0) = [0, 0, 0] and x(0) = [2, 2, −1.9] a, c or x(0) = [2, 2, −2] b, d.The red lines mark the stable periodic solution for u = 0.821 making recognition of the bifurcation more difficult, the main trend is still recognizable.
In large-dimensional nonlinear systems, it is extremely challenging to determine in advance which initial conditions will generate trajectories that pass close to the saddle-node bifurcation.Therefore, to use this technique for predicting saddle-node bifurcations, statistical quantities could be used, or general trends could be identified, and trajectories that pass far away from the saddle-node could be discarded.However, without prior knowledge of the saddle-node bifurcation, such procedures can be challenging, if not completely infeasible.
Additionally, we note that the first part of each trajectory is very hard to interpret, as it contains large quantities of all vibration modes.After the initial transient, trajectories tend to converge towards an invariant man-ifold, usually referred to as the slowest spectral submanifold [72,73] or nonlinear normal mode [74,75].Accordingly, the second part of the transient is usually much more regular, and extracting information from it is easier.Probably, implementing a modal decomposition technique [51] would solve or, at least, mitigate this problem.

Discussion and future prospects
The examples provided in Sects.3 to 6 illustrated that time series, obtained for parameter values for which there are no saddle-node bifurcations, can exhibit signs that a saddle-node bifurcation is approaching.In particular, the logarithmic decrement revealed an excellent quantity to visualize them, which can be used to predict saddle-node bifurcations in the parameter space.An important advantage of the technique is that it does not strictly require a mathematical model of the system's dynamics.Therefore, it can also be used for systems described by grey-and black-box models, whose relevance is growing in science and engineering [76,77].Besides, it can be directly implemented for experimental setup, during testing of existing devices, or for monitoring operations.
In general, time series contain a lot of information about the system dynamics, which is not easy to extract.This is proven, for example, by some recently developed system identification methods, which, based on the information of only a few time series, can extrapolate the system behavior for a significant range of parameter values [78,79].The method proposed in this study tries to extract useful information without identifying the model.A similar approach might be used to predict the occurrence of other dynamical phenomena, such as homoclinic or heteroclinic orbits, saddle-node bifurcations of equilibrium points, or Neimark-Sacker bifurcations of periodic solutions, to cite some examples.The analysis of these cases is left for future studies.
The procedure appeared to have some limitations in the case of multi-DoF systems, as illustrated in Sect.6.The main problem for large dimensional systems is that trajectories might pass very far from the region of the phase space influenced by the bifurcation to be predicted.Such trajectories present very little information about the bifurcation and are not very useful for its prediction.In these cases, obtaining several trajectories for the same set of parameter values might be necessary to reduce the probability of having only non-informative trajectories.
Another limitation of the method for multi-DoF systems is related to using the mechanical energy for the prediction.In fact, in most real systems, it is not easy to compute the mechanical energy.However, using only an approximate form of the mechanical energy may provide enough information for prediction purposes.This particular case will be investigated in future studies.
We also note that a decreasing minimum of logarithmic decrement curves does not necessarily indicate that a saddle-node bifurcation is approaching.However, from a practical perspective, it suggests that a specific region of the parameter space should be studied further to confirm whether or not a saddle-node bifurcation is occurring.
Future plans include extending the method to other bifurcations, particularly saddle-node bifurcations of equilibrium points.Besides, experimental validation of the proposed method is currently ongoing.In a wider perspective, we plan to define a test function based on a few time series, which can estimate if a system is close to a possible qualitative change of its dynamics, with vicinity interpreted both in the sense of the parameter space and of the phase space.A similar test function might be directly applied during the test phase of a device, indicating if further analysis of the system's dynamics is required.Alternatively, it could be utilized for real-time safety monitoring during operations.However, significant developments are still required for this purpose.

Conclusions
This paper presents a novel method for predicting the presence and location of saddle-node bifurcations.The method is based on measuring the logarithmic decrement, which can identify a local reduction of the energy decay of the system caused by saddle-node bifurcations.As the method relies on time series, it does not strictly require a mathematical model describing the system's dynamics, and it can be directly applied to experimental setups.The potential of the method was demonstrated by testing it on four qualitatively different examples, which confirmed its effectiveness in predicting saddle-node bifurcations.However, the study also highlighted some limitations, particularly its applicability in large dimensional systems.Nonetheless, the underlying phenomenon was present and visible in all cases.
The method's simplicity and the unnecessity of a mathematical model make it an ideal candidate for industrial applications.This possibility will be investigated in future studies.

Fig. 3
Fig. 3 Transition of a system's topology as the bifurcation parameter passes through a saddle-node bifurcation.The y-axis indicates the logarithmic decrement

Fig. 5 aFig. 6
Fig. 5 a, b, c Time series of the system in Eq. (2) for c 1 = 0.1 and various c 3 values; d, e, f corresponding phase portraits for the three different cases

Fig. 9 a
Fig. 9 a, b, c Time series of the mass-on-moving-belt system for various values of the belt velocity v; d, e, f corresponding phase portraits for the three different cases

Fig. 13 5 Fig. 14
Fig. 13 Time series of the inverted pendulum for various p values and d = 2 and τ = 0.5 Figure17adepicts the logarithmic energy decrement for various values of u.All curves present a dip that becomes increasingly pronounced as u approaches u * .Finally, the logarithmic energy decrement reaches zero as u = u * .This behavior is similar to what was observed in previous cases.However, the curves overlap and are irregular.This effect is likely related to the multi-DoF nature of the system and to the modal interaction.Although each point of the curves corresponds to a peak of the heave motion y, they are not as regular as those observed for a single DoF system.Figure17bshows a color map representing the logarithmic energy decrement for various energy and flow velocity u values.The depicted scenario exhibits a general trend similar to the previous examples.A clear reduction of the logarithmic energy decrement is visible near the saddle-node bifurcation and the periodic solutions, where it goes to zero.However, many vertical lines with different colors disturb the regularity of the color map.These lines are related to trajectories that do not pass near the saddle-node bifurcation.For topological reasons, this phenomenon cannot occur in a 2-dimensional system.Conversely, it is quite common in a multi-dimensional system, such as the present one.To further analyze this effect, we consider two trajectories obtained for the same parameter value set (u = 0.79027) with slightly different initial conditions.Figure18aand c refer to initial conditions x(0) = [2, 2, −1.9] and ẋ = 0, while Fig.18b and dto initial conditions x(0) = [2, 2, −2] and ẋ = 0.Although both trajectories converge to the trivial solution and have similar initial conditions, they exhibit very different

Fig. 16 Fig. 17
Fig.16 Time series of the system in Eq. (8) for various flow velocity values u