Stability of a nonlinear second order equation under parametric bounded noise excitation

The motivation for the following work is a structural column under dynamic axial loads with both deterministic (harmonic transmitted forces from the surrounding structure) and random (wind and/or earthquake) loading components. The bounded noise used herein is a sinusoid with an argument composed of a random (Wiener) process deviation about a mean frequency. By this approach, a noise parameter may be used to investigate the behavior through the spectrum from simple harmonic forcing, to a bounded random process with very little harmonic content. The stability of both the trivial and non-trivial stationary solutions of an axially-loaded column (which is modeled as a second order nonlinear equation) under parametric bounded noise excitation is investigated by use of Lyapunov exponents. Specifically the effect of noise magnitude, amplitude of the forcing, and damping on stability of a column is investigated. First order averaging is employed to obtain analytical approximations of the Lyapunov exponents of the trivial solution. For the non-trivial stationary solution however, the Lyapunov exponents are obtained via Monte Carlo simulation as the stability equations become analytically intractable.


Introduction
Parametrically forced second-order differential equations are of great interest in civil and mechanical engineering as they may be used to model axially loaded structures such as columns and/or pendulums. The second order equation q + 2εζωq + ω 2 1 − εμ cos νt + σW (t) + θ + εγq 2 q = 0, (1) in particular models the centre transverse displacement q of the first mode of vibration of a structural column with damping coefficient ζ, undamped natural frequency ω, and load intensity μ. The Wiener process, W (t), is used to generate a normalized bounded noise axial excitation with main frequency ν, variance σ, and initial phase θ. Finally the parameter γ captures the axial contribution to lateral stiffness of a deflected column, which is a higher order effect. The small parameter ε is inserted artificially to control the order of damping, forcing, and nonlinearity. The problem difficulty is reduced dramatically by limiting these terms to be small. Small damping is not a particularly confining assumption as most structural materials do not dissipate energy very well. The same is true for nonlinearity since for small oscillations the axial contribution to transverse stiffness is typically much smaller than bending stiffness. There is in general no reason to confine the forcing to be small, however due to the amplifying nature of resonance phenomena one can still learn a great deal about stability even with small magnitude loading.
Using variation of parameters and time scaling τ = νt, the equation of motion may be converted to two coupled first order equations of amplitude, a(τ ), and phase, φ(τ ), where κ = ω/ω 0 with ω 0 being a reference frequency. The detuning parameter Δ 0 and time scaled stochastic forcingΨ are given bȳ The deterministic case (i.e. σ = 0) is studied in detail in [1]. A similar system with viscoelastic material rather than the cubic nonlinearity is investigated for bounded noise in [2], and for wide band noise in [3]. The largest Lyapunov exponent, often simply termed Lyapunov exponent, of the amplitude, a(t), is given by It has been shown that even for stochastic systems λ is a deterministic number that represents the average growth or decay rate of the system [4]. Furthermore if λ < 0 then the system is asymptotically stable w.p.1, and if λ = 0 it is stable w.p.1. The Lyapunov exponent is only suitable for linear systems, however it may be applied to nonlinear systems by use of linearization about a solution of interest.

The trivial solution
The trivial solution, a(τ ) = 0, is of primary importance as this represents the neutral state of structural systems. The motion near the trivial solution is typically characterized by small high frequency motions about larger dominant low frequency oscillations and may therefore be approximated by averaged amplitude and phase,ā andφ respectively. This is done by the method of averaging [5]. Applying the method of averaging with κ = 1/2 to equation (2) yields equation (5), which is the first order averaged equations of motion in the resonance region. Note that for other values of κ no resonance occurs and the trivial solution is always stable.  The first order equations may be integrated directly to obtain the Lyapunov exponent of the system The expectation can be shown to be [1] E sin(2φ −Ψ) = F I (α, β) = 1 2 where α = 2Δ 0 ν σ 2 , β = − 4μω 2 νσ 2 , and I z (x) is the Bessel function of imaginary argument z. This equation may be used to investigate the affects of noise intensity, load intensity, and damping on the stability of the system. Nonlinearity does not affect stability in the vicinity of the trivial solution. Figure 2.a) shows that noise tends to improve the stability of the system since it reduces λ and narrows the region of λ > 0. The domain of λ > 0 is effectively the domain of instability of the trivial solution. Noise therefore essentially acts to disrupt resonance between the system and a periodic excitation. This has some interesting implications in the field of structural control, suggesting that artificial noise may be used to stabilize systems. Figures 2.b) and c) show that load intensity destabilizes the system and damping stabilizes the system, which is expected.

The non-trivial stationary solution
Non-trivial solutions are best characterized by the amplitude-frequency relationship, which is a plot of the amplitude of motion, a 0 , against the forcing frequency ratio, ν 2ω . In the nondeterministic case however it is necessary to use the expected value of the amplitude as the response is not a deterministic constant amplitude steady state trajectory, but rather a stationary process.
The averaged equations of motion may be used to obtain an expression for amplitudefrequency relationship of the response. Consider a stationary non-trivial stationary solution, (ā s ,φ s ) , of equation (5). By dividing the amplitude equation through byā s and taking the expected value one gets Note that the order of differentiation and integration have been reversed on the left side of both equations, which is viable since both are linear operators. Given that both of the solutions (ā s ,φ s ) are stationary processes their expection is independent of time, and therefore the time derivative of their expectation is zero. Applying definition of Δ 0 in equation (3) leads to 3.1. Special case: σ = 0 For the case σ = 0 the system is deterministic with a steady state solution, and therefore expectation operators in equation (10) may be dropped. Application of the trigonometric relationship between sine and cosine yields where the ∓ sign indicates an upper and lower solution. The lower solution is shown to be unstable in [1]. Figure 3 shows a comparison of the amplitude-frequency relationship using the simulation of the exact equations of motion in equation (2), and the analytical approximation in equation (11). There is a slight error due to the averaging procedure however the results are good, although the lower (unstable) solution cannot be simulated. Figure 4 shows the Lyapunov exponent of the trivial solution along with simulated values of the expections in equations (9) and (10). This shows that the expected value of the cosine term in equation (10) tends to zero as the noise intensity grows. The values left of the vertical dashed line are where the trivial solution has regained stability and the non-trivial solution has ceased to exist and hence they are not important. The simulation of sin term also agrees very well with the analytical result in equation (9), the small errors likely being caused by the averaging scheme. Therefore for high noise intensity the amplitude-frequency relationship in equation (10) simply becomes

Special case: high noise intensity
which implies that the upper (stable) and lower (unstable) curves converge as noise intensity increases. This coincides with the shrinking of the region of region of trivial solution instability as seen in Figure 2.

General case: Monte Carlo Simulation
The non-trivial stationary solutions for any noise intensity level may be obtained by Monte Carlo simulation. The Monte Carlo simulation comes into play in the simulation of the Wiener process for each time step. The time integration is done via Euler time stepping. The simulation needed to run for a time period long enough to allow the trajectory to converge to the non-trivial stationary solution.
Unfortunately simulation methods are only capable of obtaining stable non-trivial stationary solutions as even the smallest rounding error will cause the trajectory to diverge from an unstable solution. Figure 5 shows the amplitude-frequency relationship for changing noise intensity, load   intensity, and nonlinearity respectively. The non-trivial solution is difficult to simulate for larger frequency ratios where the trivial solution regains stability. This is because the simulations may accidentally converge to the trivial solution.

Stability
The stability of the non-trivial stationary solution is much more difficult to obtain than the trivial solution. This is because function. Choosing the reference to be the stationary solution one may write (13) The stability of the non-trivial solution is then given by the stability of the amplitude variation u.
Because the equations of variation have, as there input, the stationary solutions of the original system the stability has to be determined with a two step algorithm. The first step is to simulate a sample stationary solution. The second step is then to use that solution as the input to the equations of variation.
The general equations of variation are viable for both stable and unstable stationary nontrivial solutions. However, given that only stable solutions may be simulated, this algorithm will only work for stable solutions. Therefore if this methodology works it will always produce negative Lyapunov exponents. There is however still much to be learned from the results because it can yield relative stability (i.e. the magnitude of λ), and change in stability for changing noise, load, or damping intensity. Figure 6 shows the several Lyapunov exponent curves of the non-trivial solution for changing noise intensity, load intensity, and nonlinearity respectively. The effect of damping was difficult to study and is not included because the system converged to the non-trivial solution very slowly for small damping. Figure 6.a) shows that noise intensity tends to decrease the relative stability of the non-trivial solution in the domain of unstable trivial solution. This is the opposite affect that it has on the trivial solution. Figure 6.b) shows that while load intensity widens the band of trivial solution  instability, it does not drastically change the relative stability of the non-trivial solution at a given frequency, since all of the curves lay almost on top of each other. Finally Figure 6.c) shows that nonlinearity slightly improves the relative stability of the non-trivial solution. Figure 7 shows an overlay of typical amplitude-frequency relationships along with Lyapunov exponents of the trivial and non-trivial solutions. This shows that the peak on the right side of some of the curves is incidentally the location where the trivial solution regains stability. At this frequency ratio it was difficult to simulate a non-trivial response as most trajectories were attracted to the trivial attractor. At higher frequency ratios it became easier to simulate a nontrivial response. This is possibly due to the non-trivial response diverging from the trivial one as seen in the amplitude-frequency relationships in Figure 5. Therefore the Lyapunov exponent curves are not accurate in the vicinity of the peak on the right, however it still serves a useful role in indicating trivial solution stability. The region to the left of the first peak, is in fact the stability of the trivial solution, as it is the only solution existing in this region.

Conclusions
Lyapunov exponents are used to determine the stability of the trivial and non-trivial solutions of a nonlinear second order system subject to bounded noise parametric excitation. It is shown that damping, as expected, improves the stability of the trivial solution. Additionally noise improves the stability of the trivial solution (i.e. the preferred pre-buckled configuration), while reducing the relative stability of the non-trivial stationary solution. These results indicate that the non-deterministic components in many loading scenarios such as earthquakes might stabilize a column that would otherwise seem apt to buckle. One major drawback of the presented work is that it applies to stationary responses, which earthquakes are most certainly not. The work might, however, apply to sustained large wind loading events such as hurricanes.