Journal of Applied Mathematics and Physics
Vol.07 No.02(2019), Article ID:90522,12 pages

Random Response of a Strongly Nonlinear Oscillator with Internal and External Excitations

Claudio Floris

Department of Civil and Environmental Engineering (DICA), Politecnico di Milano, Milano, Italy

Copyright © 2019 by author(s) and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: January 9, 2019; Accepted: February 12, 2019; Published: February 15, 2019


A second order oscillator with nonlinear restoring force and nonlinear damping is considered: it is subject to both external and internal (parametric) excitations of Gaussian white noise type. The nonlinearities are chosen in such a way that the associated Fokker-Planck-Kolmogorov equation is solvable in the steady state. Different choices of some system parameters give rise to different and interesting shapes of the joint probability density function of the response, which in some cases appears to be multimodal. The problem of the determination of the power spectral density of the response is also addressed by using the true statistical linearization method.


Nonlinear Oscillator, Random Dynamics, Gaussian White Noise, Joint Probability Density Function, Power Spectral Density

1. Introduction

From the twentieth century in almost all fields of the sciences the deterministic visions have been abandoned for a probabilistic point of view. Thus, the study of the response of the dynamical system subjected to random excitations is of paramount importance. If the excitations are Gaussian white noise stochastic processes, the response of the system is a diffusive Markov vector process, whose transition probability density function (PDF) is governed by a partial differential equation that is named Fokker-Planck-Kolmogorov (FPK) equation [1] .

Unfortunately, for second and higher order dynamical systems analytical solutions of the FPK equation are known only in the final state of equilibrium or, in other words, in the stationary state of motion of the system. Moreover, if the system exhibits nonlinear damping, and if the excitations appear in the coefficients of the unknown, the so-called internal or parametric excitation, the situation is much more cumbersome. Much work has been done in this direction: see Chapter 5 of [2] and Chapter 6 of [3] and the references therein. In general, when there are a nonlinear damping and internal excitations, an analytical solution is possible only if certain relations are satisfied among the system parameters and the strengths of the exciting white noises.

This paper is aimed at studying the statistical characteristics of the response of a class of second order nonlinear systems, and highlighting the shape of response PDF by varying some system coefficients. In order to avoid using approximate methods or numerical methods, the dynamic system is chosen in a class for which the associated FPK equation admits an analytical solution [4] . When in the motion equation the damping function and/or the restoring force are nonlinear, changing some parameters may result in strange forms of the response probability function (PDF), what subtends a complicate regime of motion. To writer’s knowledge, systematic studies on this subject are lacking in literature. Thus, this paper wants to be an introductory study of this subject. This is why it has been chosen a dynamical system falling in a class for which the FPK equation is analytically solvable.

The problem of the determination of the power spectral density (PSD) of the system displacement is also addressed. While the PSD of linear systems is immediate, for nonlinear systems to writer’s knowledge only an exact solution exists [5] . Much work has been done in this field by using approximate methods [6] - [11] . Some authors assume the form of the PSD a priori, which in many instances is the PSD of a second order linear oscillator [7] [9] [10] [11] . The nonlinearities are taken into account by considering the parameters of the PSD as random quantities that depend on the amplitude of oscillation or on the system energy. The unconditioned PSD is found by integrating with respect the PDF of the amplitude or the energy.

Even in this paper the PSD of the displacement has the form of the PSD of a linear oscillator. Its parameters are found by means of the statistical linearization. As the response PDF is known in this case, it is a “true” linearization [12] , and the parameters are computed exactly. The concepts explained above are applied to some examples.

The present study has been conducted at the Technical University of Milan (Politecnico Milano), Italy.

2. Mathematical Formulation

2.1. Determination of the Response PDF

Consider the following nonlinear stochastic second order oscillator:

X ¨ ( t ) + g 0 ( X ) + g 1 ( X ) X ˙ ( t ) + g 2 ( X ) X ˙ ( t ) 2 + g 3 ( X ) X ˙ ( t ) 3 = k 1 W 1 ( t ) + k 2 X ( t ) W 2 ( t ) (1)

where the dots mean derivative with respect to time. The functions g 0 , , g 3 are deterministic ones and do not depend explicitly on time, while k 1 and k 2 are real constants. W1 and W2 are stationary Gaussian white noises with autocorrelation functions of the form E [ W j ( t ) W j ( t + τ ) ] = 2 w j δ ( τ ) ( j = 1 , 2 ) ; without loss of generality they are assumed to be uncorrelated.

The Fokker-Planck-Kolmogorov (FPK) equation associated with the dynamic system (1) is

x [ y p ( x , y , t ) ] + y { [ g 0 ( x ) + g 1 ( x ) y ( t ) + g 2 ( x ) y ( t ) 2 + g 3 ( x ) y ( t ) 3 ] p ( x , y , t ) } + ( k 1 2 w 1 + k 2 2 x 2 w 2 ) 2 p ( x , y , t ) y 2 = p ( x , y , t ) t ( y = x ˙ ) (2)

where p(x, y, t) is the joint transition probability density function of the Markov vector { X , X ˙ } . If the dynamic system is stable (to have stability, first it is sufficient that the damping is globally positive), asymptotically it tends to the equilibrium: the right-hand-side of Equation (2) becomes zero, and the reduced equation is satisfied by the equilibrium PDF. As no confusion is possible, it will equally denoted by p.

A clarification is necessary: the stochastic differential Equation (1) may interpreted in different ways (strictly speaking infinite) according to the rule of integration that is adopted, but Itô’s interpretation and Stratonovich’s one are the most important and popular. The necessity of specifying the interpretation that one adopts arises when a parametric excitation is present. However, in the present case the parametric term affects the restoring force so that Itô and Stratonovich coincide (in other words the Wong-Zakai-Stratonovich corrective terms are zero in this case [13] [14] ).

In general, Equation (2) does not admit an analytical solution even in state of equilibrium, that is with the right-hand-side null. Wang and Yasuda [4] showed that an analytical equilibrium PDF exists if the function gi have the following forms:

g 1 ( x ) = ϕ ( x ) [ c 1 2 c 2 G ( x ) ] , g 2 ( x ) = 0 , g 3 ( x ) = c 2 ϕ ( x ) (3)

where c 1 and c 2 are real constants; G(x) is a potential function, that is G ( x ) = g 0 ( x ) d x , and ϕ ( x ) = k 1 2 w 1 + k 2 2 x 2 w 2 . The above relationships impose strict conditions on the damping function, which too is made dependent on the strengths of the white noises. If the functions gi obey to the relations in (3), the equilibrium PDF is

p( x,y )=Cexp{ [ g 0 ( x )( c 1 +2 c 2 G( x ) ) ][ c 2 y 4 4 +( c 1 2 + c 2 G( x ) ) ] y 2 } (4)

where C is a normalization constant.

If the parametric excitation is absent, say k 2 = 0 , the PDF simplifies into

p ( x , y ) = C exp [ 1 k 1 2 w 1 ( c 1 Λ + c 2 Λ 2 ) ] (5)

where Λ = y 2 / 2 + G ( x ) is the mechanical energy of the oscillator. It is noted that Equation (5) is also the PDF of the energy, if this is considered an independent variable.

2.2. Spectral Characteristics of the Response

As advanced in the Introduction, approximately it is assumed that the response power spectral density (PSD) has the same form as that of a linear oscillator, that is

S X X ( ω ) = k 1 2 w 1 ( ω e 2 ω 2 ) 2 β e 2 ω 2 (6)

In Equation (6) the parameters of the PSD ω e , β e will be computed by means of the statistical linearization method [15] . As this method is not suitable when there is a parametric excitation, the case of only external excitation is considered, whose strength appears in the numerator of Equation (6).

The nonlinear system (1) is replaced by the following linear system equivalent to (1) with k 2 = 0 in some statistical sense:

X ¨ ( t ) + β e X ˙ ( t ) + ω e 2 X ( t ) = k 1 W 1 ( t ) . (7)

Using Equation (7) instead of (1) causes the following error:

E = g 0 ( X ) + 2 k 1 2 c 2 w 1 G ( X ) X ˙ + k 1 2 c 2 w 1 + β X ˙ β e X ˙ ω e 2 X (8)

where β = k 1 2 c 1 w 1 . Minimizing the error E in mean square, one obtains

β e = β + E [ z X ˙ ] E [ X ˙ 2 ] , ω e 2 = E [ z X ] E [ X 2 ] (9)

where z = g 0 ( X ) + 2 k 1 2 c 2 w 1 G 0 ( X ) X ˙ + k 1 2 c 2 w 1 . The question arises as the statistical averages in (9) are to be evaluated. In [7] and [11] they are considered amplitude dependent and the dependence is eliminated by averaging Equation (6) with respect to the PDF of the amplitude. In the present case, the PDF of the response is known exactly so that the averages are computed with respect to it (the computations are numerical, but even in the papers cited among the references the computations are so). When the linearization is done by using the exact PDF, it is called true linearization [12] .

3. Applications

3.1. Dynamic Systems with Internal and External Excitations

The function g 0 ( X ) represents the restoring force of the dynamic system: it has been chosen

g 0 ( x ) = d G ( x ) d x = a x + b x 3 , G ( x ) = 1 2 a x 2 + 1 4 b x 4 (10)

The constant b must be positive, otherwise the system diverges, while a my be negative or positive. If it is positive, the potential function G has only a minimum in the origin. When it is negative, G has an unstable maximum in the origin and two minima in ± a / b (Figure 1 left), which causes complicated dynamics. In the analyses a assumes the values 1 , 0 , + 1 .

The function F d = ( k 1 2 w 1 + k 2 2 w 2 x 2 ) ( c 1 + 2 c 2 G ( x ) ) x ˙ + c 2 ( k 1 2 w 1 + k 2 2 w 2 x 2 ) x ˙ 3 is the damping function of the system. If c 2 < 1 , there is an interval in which the damping is negative (Figure 1 right), that is energy is introduced in the system. In the analyses the values c 2 = 2 and c 2 = 1 have been considered.

All the computations and the plots have been performed by using the software MAPLE ruled on a work-station. The computational charges have been low in every case: at the most few seconds.

In Figure 2 there are the joint PDF of the case a = 1 , c 2 = 1 and its sections x ˙ = y = 0 , x = 0 (Equation (4)). The other parameters take the values: b = 0.5 , k 1 = 1 , k 2 = 0.5 , w 1 = w 2 = 1 , c 1 = 1 . The joint PDF resembles a sand cone like a Gaussian bimodal PDF. Surely, the tails of this PDF are more pronounced than in a Gaussian, but this type of plot cannot show this fact. Due to the regularity of the PDF, no other comments are necessary.

Figure 3 reports the same graphs as in previous figure for the case a = 0 (the other parameters are unchanged). It is evident that the PDF is rather flat in x-direction, which is confirmed by the section x ˙ = y = 0 , so that it does not resemble a Gaussian PDF. The kurtosis E [ X 4 ] / σ X 4 confirms that this case is well far from the Gaussian one: its value is 2.066 against 3 of a Gaussian PDF (it is evaluated numerically). On the contrary, the section x = 0 is very regular and Gaussian like.

Things are very different when the parameter a in Equation (10) becomes negative: a = −1 has been chosen (the other parameters are unchanged). The potential function is now double-well (see Figure 1 left). If there is sufficient energy, the displacement continuously moves from a well to another. This fact is reflected by the shape of the PDF (Figure 4): it has two modes, and the sections

Figure 1. Potential function: a < 0 , left. Dissipation function with c 2 = 2 ( y = x ˙ ) , right.

Figure 2. Joint PDF of the system with external and internal excitation (Equation (4)), a = 1. Top: three-dimensional plot. Bottom: section x ˙ = y = 0 (left), section x = 0 (right).

Figure 3. Joint PDF of the system with external and internal excitation (Equation (4)), a = 0. Top: three-dimensional plot. Bottom: section x ˙ = y = 0 (left), section x = 0 (right).

Figure 4. Joint PDF of the system with external and internal excitation, a = 1 (Equation (4)). Top: three-dimensional plot. Bottom: section x ˙ = y = 0 (left), section x = 0 (right).

in x-direction are always bimodal, while the sections in orthogonal direction are always unimodal Gaussian like. It is not new that the PDF of the displacement is bimodal in the case of double-well potential: as there are many and many papers that report this fact, we renounce to cite them. However, it is emphasized that in the present case in x-direction the PDF is rather flat and the saddle is not deep: see Figure 4, bottom, left.

As advanced at the beginning of this section, the damping function is negative over a certain interval when c 2 < 1 : it has been chosen c 2 = 2 in the next two examples, while a may be −1 or +1 with the other parameters unchanged.

The negative damping even in a short interval changes the things dramatically. The overall aspect of the PDF is crater-like, but a cone is inner as it is revealed by the section x ˙ = 0 , which has three modes (Figure 5). To writer’s knowledge, the examples of crater-like PDF’s are not numerous in literature (see the example on page 185 of [2] and [16] [17] ), but probably a crater-like PDF with a maximum in its interior has not yet presented under stochastic excitations only. Three-modal PDF’s are very rare. In [16] the system is acted by a harmonic excitation and by a Gaussian with noise having an intensity small with respect to the intensity of the harmonic excitation; moreover, the nonlinear term too is assumed to be small. The sections x = constant may exhibit two peaks as in right plot at the bottom of Figure 5 or only one peak (for brevity’s sake the section

Figure 5. Joint PDF of the system with external and internal excitation, a = 1 , c 2 = 2 (Equation (4)). Top: three-dimensional plot. Bottom: section x ˙ = y = 0 (left), section x = 0 (right).

x = 1 is not shown). Another characteristic that has not been plotted is that the level curves (p = constant) tend to be angular nearly rectangular.

In next application of this section the parameters are a = 1 , c 2 = 2 . The three-dimensional plot of the PDF and its sections are in Figure 6. Still the PDF is crater-like but with a simpler and more regular aspect. In both the sections x ˙ = 0 and x = 0 the saddles are rather deep. However, with a sharper examination of the three-dimensional plot it can be concluded that in both x and x ˙ there are unimodal sections.

3.2. Dynamic Systems with Internal Excitation Only

In the case of external excitation only the constant k 2 is zero: for clarity’s sake we report the motion equation of the system, which is

X ¨ ( t ) + g 0 ( X ) + k 1 2 w 1 [ c 1 + 2 c 2 G ( X ) ] X ˙ ( t ) + k 1 2 w 1 c 2 X ˙ ( t ) 3 = k 1 W 1 ( t ) (11)

The joint PDF of X and X ˙ is given by Equation (5). The response PSD is computed by using the method of Section 2.2. In the analyses a is worth −1, while c 2 takes the values 1, −2. The other parameters are: b = 0.5 , c 1 = 1 , k 1 = 1 , w 1 = 1 .

In Figure 7 there are the results of the case a = 1 , c 2 = 1 . The joint PDF is bimodal in x-direction as the potential is double-well, while in the direction of the velocity is unimodal. It is recalled that the power spectral density of the

Figure 6. Joint PDF of the system with external and internal excitation, a = 1 , c 2 = 2 (Equation (4)). Top: three-dimensional plot. Bottom: section x ˙ = y = 0 (left), section x = 0 (right).

Figure 7. Joint PDF (Equation (5)) of the system with internal excitation of Equation (11), a = 1 , c 2 = 1 (left). Power spectral density of the system according to the true linearization method (right).

displacement X(t) is defined as

S X X ( ω ) = 1 2 π ω + ω R X X ( τ ) e i ω τ d τ (12)

where R X X ( τ ) is the autocorrelation function of X(t), and i = 1 (using the method of Sec. 2.2, it is not necessary to know R X X ( τ ) . Thus, S X X ( τ ) extends on the whole real axis, and it is symmetric with respect to the vertical axis in the origin. It is expressed by Equation (6), in which the parameters are computed by means of the formulae in Equation (9): ω e = 0.9442 , β e = 3.6674 . Being the linearized frequency ω e less than one, the PSD has a maximum only, which is in the origin.

Figure 8 reports the results of the case a = 1 , c 2 = 1 . For this value of c 2 there is an interval in which the damping is negative, that causes a crater-like

Figure 8. Joint PDF (Equation (5)) of the system with internal excitation of Equation (11), a = 1 , c 2 = 2 (left). Power spectral density of the system according to the true linearization method (right).

1It is recalled that the nonlinear system method allows replacing a dynamic system for which the FPK equation is no solvable by means of one for which the FPK is analytically solvable. However, the computations are generally cumbersome.

PDF as in Figure 5 and Figure 6. The power spectral density has two symmetric peaks that are apart from the origin. On the other hand the value of the PSD in the origin is not negligible so that the displacement has energy in a long interval of frequency. It might be contested that the power spectral densities of nonlinear systems have more peaks. The method adopted herein may be viewed as a first order method so that it is not capable of detecting secondary. However, if secondary peaks are detected [10] [11] , these are not important. In this case, the linearization parameters are: ω e = 5.2682 , β e = 3.8453 .

4. Concluding Remarks

In this paper, the statistical characteristics of strongly nonlinear second order oscillators are studied by examining the joint probability density function of the response X , X ˙ . It has been chosen a class of oscillators for which the Fokker-Planck equation is exactly solvable in the equilibrium regime. In the general case, the excitation is formed by an external stationary Gaussian white noise and by an internal (parametric) Gaussian white noise proportional to the displacement. But, in order to compute the power spectral density function of the displacement the internal excitation is absent as the method of analysis, statistical linearization, is defective in the presence of internal excitations.

Both the damping function and the restoring force are of polynomial form, being the polynomials defined by some real constants. For the Fokker-Planck equation having an analytical solution, strict relationships must be satisfied, which link the strengths of the white noises with the system parameters. This fact was already known [2] [3] , and constitutes a serious limitation in using the Fokker-Planck equation method as in real systems the necessary relationships are rarely met1. The rationale of choosing such a type of system is to avoid using approximate methods or numerical ones. In fact, the paper is aimed to study as the response probability density function varies when the system parameters are varied. The applications show that the response statistics vary substantially by varying the constants that define the damping function and the restoring force of the system. Crater-like multimodal PDFs are detected that ceil complicated dynamics. To writer’s knowledge, such a type of studies are lacking in literature. Other systematic analyses are surely necessary.

The dynamic characteristics of a system are difficult to be estimated, and in practice they are adapted to a linear model, whose response might be very far from the actual response, if in reality the oscillator is nonlinear. The lesson that can be drawn is that one must be very cautious in choosing a linear model.

The determination of the power spectral density of a nonlinear system is still an open problem, for which the studies are not numerous [5] - [11] . In most cases the form of the spectral density is chosen a priori as that of a linear oscillator, but with random parameters. This is done in this paper too. As the FPK equation is exactly solvable in the cases considered here, the random parameters of the spectral density are obtained by averaging with respect to the exact probability density function, which results in simpler calculations with respect to the methods of [6] [10] [11] . However, the present approach does not allow detecting secondary peaks in the spectral density, and it is not applicable when the response PDF is not known analytically.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

Cite this paper

Floris, C. (2019) Random Response of a Strongly Nonlinear Oscillator with Internal and External Excitations. Journal of Applied Mathematics and Physics, 7, 331-342.


  1. 1. Risken, H. (1996) The Fokker-Planck Equation. Methods of Solution and Applications. Springer, Berlin.

  2. 2. Lin, Y.K. and Cai, G.-Q. (1995) Probabilistic Structural Dynamics: Advanced Theory and Applications. McGraw-Hill, Inc., New York.

  3. 3. Cai, G.-Q. and Zhu, W.-Q. (2017) Elements of Stochastic Dynamics. World Scientific, Singapore.

  4. 4. Wang, R. and Yasuda, Z. (1998) Exact Stationary Solution of Six Classes of Nonlinear Stochastic Systems under Stochastic Parametric and External Excitations. ASCE Journal of Engineering Mechanics, 142, 18-23.

  5. 5. Dimentberg, M.F. (1995) Spectral Density of a Non-Linear Single-Degree-of-Freedom System’s Response to a White Noise Random Excitation: A Unique Case of an Exact Solution. International Journal of Non-Linear Mechanics, 30, 673-676.

  6. 6. Dykman, M.I. and Krivloglaz, M.A. (1985) Spectral Distribution of a Nonlinear Oscillator Performing Brownian Motion in a Double-Well Potential. Physica A, 133, 54-73.

  7. 7. Miles, R.N. (1989) An Approximate Solution for the Spectral Response of Duffing’s Oscillator with Random Input. Journal of Sound and Vibration, 132, 43-49.

  8. 8. Valéry Roy, R. and Spanos, P.D. (1993) Power Spectral Density of Nonlinear System Response: The Recursion Method. ASME Journal of Applied Mechanics, 60, 358-365.

  9. 9. Bouc, R. (1994) The Power Spectral Density of Response for a Strongly Non-Linear Random Oscillator. Journal of Sound and Vibration, 175, 317-331.

  10. 10. Krenk, S., Lin, Y.K. and Rüdinger, F. (2002) Effective System Properties and Spectral Density in Random Vibration with Parametric Excitation. ASME Journal of Applied Mechanics, 69, 161-170.

  11. 11. Spanos, P.D., Kougioumtzoglou, I.A. and Soize, C. (2011) On the Determination of the Power Spectrum of Randomly Excited Oscillators via Stochastic Averaging: An Alternative Perspective. Probabilistic Engineering Mechanics, 26, 10-15.

  12. 12. Kozin, F. (1987) The Method of Statistical Linearization for Non-Linear Stochastic Vibrations. In: Ziegler, F. and Schueller, G.I., Eds., Nonlinear Stochastic Dynamics Engineering Systems, Springer, Berlin, 45-56.

  13. 13. Wong, E. and Zakai, M. (1965) On the Relation between Ordinary and Stochastic Differential Equations. International Journal of Engineering Science, 3, 213-229.

  14. 14. Stratonovich, R.L. (1966) A New Representation for Stochastic Integrals and Equations. SIAM Journal of Control, 4, 362-371.

  15. 15. Roberts, J.B. and Spanos, P.D. (1990) Random Vibration and Statistical Linearization. John Wiley & Sons, Chichester, UK.

  16. 16. Von Wagner, U. (2002) On Double Crater-Like Probability Density Functions of a Duffing Oscillator Subjected to Harmonic and Stochastic Excitation. Nonlinear Dynamics, 28, 343-355.

  17. 17. Feng, G., Xu, W., Rong, H. and Wang, R. (2009) Stochastic Responses of Duffing-Van der Pol Vibro-Impact System under Additive and Multiplicative Random Excitations. International Journal of Non-Linear Mechanics, 44, 51-57.