Prediction of Extreme Response of Nonlinear Oscillators Subjected to Random Loading Using the Path Integral Solution Technique

This paper studies the applicability of the path integral solution technique for estimating extreme response of nonlinear dynamic oscillators whose equations of motion can be modelled by the use of Itô stochastic differential equations. The state vector process associated with such a model is generally a diffusion process, and the probability density function of the state vector thus satisfies the Fokker-Planck-Kolmogorov equation. It is shown that the path integral solution technique combined with an appropriate numerical scheme constitutes a powerful method for solving the Fokker-Planck Kolmogorov equation with natural boundary conditions. With the calculated probability density function of the state vector in hand, one can proceed to calculate the required quantities for estimating extreme response. The proposed method distinguishes itself by remarkably high accuracy and numerical robustness. These features are highlighted by application to example studies of nonlinear oscillators excited by white noise.


Introduction
An important element in the safety assessment of many engineering systems, is the task of estimating the probability of extreme events that may Jeopardize the structure in some specified sense.Very often, this problem can be formulated as finding the probability that some time varying random quantity docs not exceed a specified capacity level during a given time period.Stated this way, the problem typically reduces to a study of the extreme values of a .stochasticprocess originating as the response of a system subjected to some stochastic loading process.
In this paper the focus will be on the problem of estimating the extreme response of nonlinear dynamic systems subjected to random forcing pnx;esses.In recent years the methods of time domain Monte Carlo simulations, see, e.g., Refs.fl-S], have received consid-erable attention as a tool for estimating response statistics.These methcxls are versatile and attractive in the sense that nonlineariiies can be easily dealt with.TTie main drawback at present is the large CPU times needed for accurate prediction of extreme responses.Even if this issue seems to become less of an obstacle every year, portending perhaps that such methods may dominate practical estimation of response statistics of nonlinear systems in the not too distant future, it will still be desirabk to ha\'e available alternative methods of calculating the response statistics, both simplified and more elab<')rate.Here we shall explore a method ba.scd on the theory of Markov diffusion processes.The justification for using this theory is related to the fact that the response of nonlinear dynamic systems to broad band random excitation can very often be accurately dc-scribed by applying the theory of multidimensional Markov priKCSses.By this, the extensive theory of Markov diffusion processes can be brtxight to bear on these problems.In particular, it can be shovi^n that the probability law of response quantities can be derived by solving a partial differential equation, viz., the Fokker-Planck (-Kolmogorov) (FPK) equation, see Refs.[6.7].In most cases of practical interest, this equation has to be solved numerically.
In the next section we shall describe a method for solving the FPK equation that is ba.scd on a formal solution of the same equation.This solution is obtained by invoking the fact that a Markov diffusion process locally looks like a Brownian motion.By using the Markov property, the global solution can dien be constructed by linking the local solutions, which are known explicitly.The obtained solution is generally known as a path integral solution (PIS).The reader is referred to Ref. [7] for a further discusskm.One of the first efforts to cxptoit the PIS method explicitly in developing numerical solution algorithms is described in Ref. [81.Subsequently, other authors have also used the PIS approach to solve various random vibration problems, cf.Rcfs.[9- ! 3].
Before embarking on a description of the PIS method, it is expedient to briefly show how the obtained solutions are used in an extreme value analysis.Assuming that the response quantity of interest is a scalar (real) stationary stochastic process, Z(t) say, the PIS method typically provides a numerical estimate of the joint probability density function (PDF)/ii(».«) of Z(r) and Z(f)=dZ(r)/ dr.It is now assumed diat the mean level upcrossing rate vj(') of Z(f) can be calculated from Rice's formula a.s follows ^{z)-\yfzz(z,y)dy. (1) Adopting the assumption that upcrossing of high levels are statistically independent events, which leads to Poisson distributed crossings, it follows that an asymptotic approximation of the probability distribution function of the extreme value of the process Z(t) during a time T, denoted by M(T) (=sup{Z(/);0:£/<r}). is given by Prob{Ma)-^z}-cxpi-i^Xz)T) (7"^^). ( The accuracy of Eq. ( 2) depends to a large extent on the effective bandwidth of the response process Z(f).Decreasing bandwidth leads eventually to a significant clumping effect of large response peaks, invalidating the assumptkin of statistically independent upcrossing of high levels.Methods that aim at correcting for this effect have been proposed for Gaussian (Refs.[14,151) and non-Gaus.sian(Ref.[ 16]) processes.However, thus poini.will not be pursued any further here.We shall assume that Eq, (2) provides an acceptable approximatkm.Hence, the central parameter to be determined is the upcrossing frequency i^K*), which Ls easily calculated once the joint PDF/zz<*.') of Z{r) andZ(f)-dZ(f )/df has been made available.In the next section it is shown how this PDF can be calculated for the resEK>nsc of a wide range of nonlinear oscillators subjected to white noise or filtered white noise loading.
(3) has the final form to be used subsequently.
It is demonstrated in Ref.
[6], it can be proved that the TPD p(x,t \x',i') (f^r'aO) is the solution of a partial differential equation of the form O denote appropriate zero-matrices and G(*) denotes an (n-r)X(rt-r)-matrix function with elements g,j('), jj=r+l,...,n.G(*) will be called the reduced diffusi<in matrix.Equation ( 6) can now be rewritten as where Let/(jc.r) denote the PDF of the random vector ^(r).U f(x,!')=w{x) for some initial PDF w(x), then it is recognized from Eq. ( 6) and the relation lhat/{je,/) itself is a solution of Eq. ( 6) satisfying the initial condition/(X,;')=K'(X).
In this paper wc shall be interested primarily in stationary .solutions/,(j) to [iq.(6), that is /,(xHim/(x.f)=Iim/?(Jr,/\x',i') (9) provided they exist.Even when both limits exist, it is clear that lim/(x,/) provides the faster convergence when the initial condition/(x,f')*=/,(x).This comment is relevant to the numerical implementation of the PI.S method, and will be discussed below.11) obtain the PIS appropriate for the dynamic systems studied in this paper, it is necessary to be more specific on the structure of the matrix function Q{').In particular, it will be assumed that the first r rows of Q(') arc zero, that is q.ji')^ for '=1 '■; j=U..J>i (r<n) (10) and that ^,/*)#0 for at least one j for every iW+1 n.This implies that the diffusion matrix G(*) assumes the form Proceeding in a manner similar to the derivations in Ref. [7], it can be shown that the TPD for small values of T(-/-r') is given by the following expression, which is correct up to terms of order T^ p{x,t+T where I GI denotes Ihc determinant of the reduced diffusion matrix G, assumed to be positive definite.This implies that 161 > 0. [6" '],j denotes the element in position ij of the inverse matrix of G.As shown in Ref. [7], the expression given by Eq. ( 13) is not unique, but seems to be well suited for tnir purpose.
Having obtained an explicit expression for the TPD for a short time step, one can now invoke the Markov property.This allows a TPD over a time interval of arbitrary length to be expressed in terms of a product of short-time TPDs.% dividing a given time interval (r',/) into A' small time intervals of length T^t-t'VN, it is found that (f,=?'+/T, /=f \ x-x''^\ f'=f".x'=x"") p(xAx\t')=\-{f\p{x^jt,\x^''\t^^W (14) Similarly, with an initial PDF/(x'.f')=w(x').the PDF f{x4) will be given by /(x .t)= f -f n P (x''\lj IX*' '',/, , )H'{x"'Odi'"'..nlx''' ''.(15) Hence, by combining Eq. ( 13) with Eqs. ( 14) or (15), a formal (approximate) solution of the FPK equation can be written.Equations ( 14) and (15), which arc often referred to as PIS, constitute the core of the numerical solution procedure to be described subsequently.It is realized chat a numerical solution according to this method, automatically provides the evolution in time of the (conditional) PDF of the Markov process Xit) from given start conditions in terms of an initial density fix' ,t')=w(x'), including the degenerate case /(jr*,r')-S(jr'-Jro), for some starting point Xo.It is also worth noting how the PIS relates to the physics of the dynamic model, which is expressed through the coefficients ntji*) and £?,;(•), cf.Eq. (3).The evolution in time of the PDF as expressed by the PIS, is seen to be directly determined by these coefficients in an explicit manner.This fact is a very important advantage of the PIS method, and reveals its fundamental physical significance.

Numerical Implementation
In the numerical implementation, the PIS is obtained by an iteration process based on the Chapman-Kolmogorov equation expressed as The discretization of state space for the numerical solution makes it appropriate to employ an interpolation and smoothing procedure to increase the numerical efficiency.It was found that application of cubic B-splines, as detailed in Ref. [17), offered the desired accuracy and smoothness for the type of problems considered in this paper.This procedure was used as follows.At each time step tj-,-itj,p(x^~'\t,-i \x\t') is represented as a cubic B-spline series in the following manner

Examples
The accuracy and power of the developed PIS prt^edurc will be illustrated by application to specific case studies taken from two classes of dynamic models.Both models are described by Eq. ( 3) with n-2 and m-3.This implies a two-dimensional state space vector X={X,J(2f'-{ZJZf.Further, 171 (•) and Q(*) are such that '"[(A'lA)-^: and q,j{')-0 for_/=l,2,3.Assuming sufficient restrictions on m{») and Q('), cf.Refs.[6,7], X{t) becomes a Markov diffusion process.Invoking Eq. ( 13), it can be shown that, up to correction terms of order T^, the associated TPD assumes the form By combining Eqs. ( 25) and (26), aixi applying (he solution technique described in the previous section, the TPD p(x,t I jc',f') for large t-r' can be calculated.By this, the time evolution of the system when it starts from rest, for example, can be studied.The stationary PDF is obtained in the limit as f-f'->^.For application of the PIS method to other problems involving both two-and three-dimensional slate space vectors, the reader may consult Reft.111-13,18.19],

Example 1-The Caughey Oscillator
There is a class of dynamic models for which there exist an analytical solution for the stationary joint PDF of X.A member of this class m^ be called a Caughey oscillator, Ref.
[20], The generic equation of motion for this oscillator can be written as N(t) denotes a stationary, zero-mean Gaussian white noise satisfying E[N{t) A'(r+T)]-B(T), where &(•) denotes Dirac's delta function, /'is a positive con.stanl and g(E) is a function of the total energy E=E(Z^) given as follows where For this example m2(z4)~ zs[E(.z,z)]-h(z),and we may choose qi\=qi2='W,-W2=0.qii-F and dW,{,t)-N{t)dt.The stationary, joint PDF, denoted by p^{*), is then determined by the relation, cf, Refs.
For the illustration purposes in this paper, we hjtve chosen the following special case of Eq. ( 28) with parameters ^, e, and A.
The stationary PDF only depends on the parameters e and A, and the numerical solution for the following set of parameter values has been calculated (c,A)=(0.0)(Gaus.sianresponse), (0, 0.2) (Duffing oscillator) and (0.5, 0.1).The calculations were carried out with the same number of grid points on both axes in state space, aviz-, 45.Since the resulting PDFs are actually independent of ^, the value §=0.\ was chosen for the Gaussian and Duffing ca.ses, while ^-0.5 was adopted for the last case.The time increments used were T-0.0025 S, 0.001 s, and 0.02 s, respectively.The total CPU time on a DEC station 3100' was about 5 minutes for each case.In Figs.. 1 and 2 are shown the marginal PDFs of the displacement response for the three case studies considered, together with the corresponding analytical solutions, in Fig. 3 are given the corresponding analytical and numerical results for the mean upcrossing rate.It is seen that in all three cases the agreement between the numerical PIS and the analytical solution is very gotxl over the whole range of probability levels given.In fact, the accuracy can be retained down to much lower probability levels (=10 '") at a mtxJcratc increase in computer time.
This model was studied by Dimentberg [22], who showed that when w,n=^en In this example, the resptmse statistics of a nonlinear oscillator subjected to both external and parametric random excitation will be illustrated by applying the methodology of the paper to two specific case studies.
The equation of motion of the oscillator is the following Z+2^[l+;V,(/)]Z+7[Z'+-TlZ+ta?ri+^2(/)]Z-iV3(r).By this, we have the opportunity to test the accuracy of the PIS method for this kind of dynamic model.The results of two particular cases will be presented.Case 1: Here the following parameter values were used.^.1, T^.l, wi,-1.0,/T=2.5, /7=0.l,A?-0.3.For the numerical calculations a grid size of 49X49 points and a time increment 7M).01 s was used.The total CPU time on a DEC 3100 work station was 3 min for the PIS calculation.The results for the analytical and numerical solutions are given in Figs.4-6

Conclusions
A numerical method for estimating the extreme response of nonlinear oscillators excited by white noise, or filtered white noise, has been described.The example calculations presented show that the method gives very accurate estimates of the required joint PDF.In fact, for every example having analytical solution on which the method has been tested, complete agreement has been found with proper choice of grid size and time increment in the numerical solution procedure.In the present paper, of course, only a few cases can be given.Experience with the method indicates that two-dimensional problems can be solved routinely with high accuracy requiring a few minutes CPU time on a work station (DEC .station3100).The solution of three-dimensk)nal problems requires more care in the sense that computer capacity starts to become an issue of importance.In such cases the CPU time easily runs into htxirs.