Numeric Solution of the Fokker-Planck-Kolmogorov Equation

The solution of an n-dimensional stochastic differential equation driven by Gaussian white noises is a Markov vector. In this way, the transition joint probability density function (JPDF) of this vector is given by a deterministic parabolic partial differential equation, the so-called Fokker-Planck-Kolmogorov (FPK) equation. There exist few exact solutions of this equation so that the analyst must resort to approximate or numerical procedures. The finite element method (FE) is among the latter, and is reviewed in this paper. Suitable computer codes are written for the two fundamental versions of the FE method, the Bubnov-Galerkin and the Petrov-Galerkin method. In order to reduce the computational effort, which is to reduce the number of nodal points, the following refinements to the method are proposed: 1) exponential (Gaussian) weighting functions different from the shape functions are tested; 2) quadratic and cubic splines are used to interpolate the nodal values that are known in a limited number of points. In the applications, the transient state is studied for first order systems only, while for second order systems, the steady-state JPDF is determined, and it is compared with exact solutions or with simulative solutions: a very good agreement is found.


Introduction
It is widely recognized that many types of agencies acting on engineering structures and equipments have random characteristics.This is the case of earthquakes, wind load, wave load, electric current in a circuit, and so on.Thus, these agencies are to be described by means of the theory of stochastic processes.Moreover, in some cases, the structural properties themselves are uncertain, and must be considered as stochastic processes, which gives raise to a multiplicative or parametric excitation.In this way, the differential equations that govern the response become stochastic differential equations.
If the differential equations that govern the response are nonlinear or the excitation is parametric or both, the random response analysis is not an easy task, and mathematically exact solutions for the statistical characterization are not available in many cases.However, if the excitations are Gaussian white noise (WN) random processes, the system response is a vector Markov process, too: see [1][2][3].
If the initial state of a Markov process is known, then the joint probability density function (JPDF) of the states characterizes completely the stochastic process.The JPDF is the solution of a parabolic partial differential equation (PDE), the so-called Fokker-Planck-Kolmogorov (FPK) equation.This equation depends on time and on the actual values of the system states.However, if the system is damped and stable, a stationary state exists, which is characterized by the stationary JPDF p(x), which does not depend on time.The FPK equation loses the dependence on time, and is called reduced FPK equation.A wide treatment of it can be found in [3,4].The severe limitations of the FPK equation approach resides in that exact analytical solutions are available in a restricted number of cases, and mostly in the steady state.In the transient state, analytical solutions were obtained for scalar systems only: the reader is referred to [4].Along the years, much theoretical work has been done in order to solve the reduced FPK equation: .Notwithstanding the considerable progress in this matter, two serious flaws remain: 1) in the case of multiplicative excitations, for the FPK equation to be solvable, restrictive relationships among the system parameters and the spectral densities of the excitations must be satisfied, and such requirements are rarely encountered in practical cases; 2) in several cases, no analytical solutions have been found, among the case of oscillators with a generic nonlinear damping mechanism.
For these reasons, the sixties numerical schemes of solution have been developed parallelly to the theoretical studies.These methods are: weighted residual method, eigenfunction expansion, finite differences, variational principles, and finite element (FE) method.
In the weighted residual method [28][29][30][31][32][33], the functional form of the JPDF is chosen a priori: a mixture of Gaussian functions or the exponential of a polynomial in the state variables is the most usual choice in both unknown coefficients' appearance.The approximate form of the JPDF cannot satisfy the FPK equation exactly so that a residual error arises.It is imposed that the projection of this into a set of weighting functions vanishes.In this way, a system of equations in the unknown coefficients is obtained, and it must be solved.Muscolino, Ricciardi and Vasta's method [29] is notable as it gives raise to an expression of the JPDF in which the coefficients are linear functions of the response statistical moments.
In the eigenfunction expansion method [34][35][36][37][38], the non-stationary JPDF of the system states is expressed as a truncated series of orthogonal functions, which are functions that form a complete orthonormal basis in a Hilbert space [39].If the eigenfunctions of the FPK equation were available in all cases, this representation would be exact, but in general they are not known, and their computation even in an approximate numerical form is often cumbersome.Thus, the different authors adopt different strategies: in [35,38], the orthogonal functions are chosen a priori, and in order to determine the coefficients of the series, the weighted residual method is used.In [34,36,37], a perturbative approach is used.
Roberts use the finite difference method [40] to solve the FPK equation for the one-dimensional PDF of the energy envelope.In other words, there is only a spatial variable beyond the time variable, which is a limitation.
In [41], two types of variational approaches are presented, both of which are aimed at finding approximate values of the non-zero eigenvalues of the FPK operator (the first eigenvalue is zero, and it corresponds to the steady state PDF).The former one constructs an approximate Hermitian operator, and the other is based on a perturbation expansion.The applications are confined to the steady state.In [42], the variational approach gives raise to an iterative procedure.The applications regard the transient state of scalar systems.
The approximate solution of PDE's by means of the finite element method (FE) is a classic topics in numeric mathematics (e.g.see [43]).Probably, Bergman and Heinrich [44][45][46] were the first who applied this method in the field of stochastic dynamics.In the above cited references, they did not analyze the FPK equation but the Pontriagin-Vitt equation in the moments of the first passage time of a level x from the displacement X of a second order oscillator (as regards the Pontriagin-Vitt equation see [3], Chap.8); however, the principles are the same.Then, the FE method was applied to the FPK equation: [47][48][49][50][51][52][53][54].Differently from the weighted residual method in which the transition JPDF of the state vector has the same approximate form in the whole state space, in the FE method, the state space is discretized into elements inside which the JPDF varies according to the shape or trial functions that have been selected previously, and depending on the nodal values (see later).In [47][48][49][50][51][52][53][54], linear shape functions are used, which insure C 0 continuity only.Vice versa, as the FPK equation Equations ( 6) and ( 7) contains second derivatives, it would require C 1 continuity.The problem is overcomed by combining the FE method with the weighted residual method: the residual error is made orthogonal to a weight function in such a way that the coefficients of a linear system in the nodal values are computed.Fundamentally, there are two versions of the FE method for solving the FPK equation: the Bubnov-Galerkin and the Petrov-Galerkin method.In the former, the weighting functions are equal to the trial functions.Vice versa, in the latter, weighting and trial functions are different, in general the weighting functions being non linear.This version is believed more suitable for convective problems [44,45,49,52].However, in the steady state, that is when solving the reduced FPK equation, there is no diffusion of probability.With the exception of [49,52,53] in the above mentioned references, the applications regard the steady state because of the charge for computation.Some improvements of the method were proposed: in [51], an adaptive grid generation is used; in [53], a multi-scale FE method is fitted to the present problem; in [54], the solution is improved by means of a local hp-refinement of the mesh.
The present paper is organized as follows: Sec. 2 is devoted to the FPK equation.Sec. 3 treats the FE method for solving the FPK equation: first, its theoretical derivation is reviewed, and some topics are discussed.In detail, it is shown that: 1) to interpolate the nodal, values with splines of degree higher than those of the shape functions may improve the solution substantially; 2) in some cases the use of weighting functions different from the shape functions yields a notable computing time saving.The last two sections are devoted to the applications and conclusions respectively.In order to contain the computing time, the transient state is studied for first order systems only, while for second order systems the steady-state JPDF is determined.

The Fokker-Planck-Kolmogorov Equation
Let     g X t an arbitrary function of a stochastic process X(t), which at its turn is solution to the stochastic differential equation where B(t) is a unit Wiener process, is the drift term, and is the diffusion term.As the latter depends on X(t), that is the excitation is multiplicative, it is intended that includes the so-called Wong-Zakai-Stratonovich corrective term [55,56], which reads as where is the originary drift term before correction.
Now, we apply Itô's differential rule [57][58][59] to the first derivative of the stochastic expectation of g: Expressing the expectations with respect the conditional PDF  0 0 , , Integrating Equation (4) by parts, and taking into account that the PDF p and its first derivative tend to zero as x tends to infinity; we obtain As g(x) is quite arbitrary, it may be assumed   1 g x  , which yields: Equation ( 6) is the Fokker-Planck-Kolmogorov equation (or forward Kolmogorov equation) in the unknown transition JPDF of the scalar process X(t In Equation ( 7) the summation rule with respect to repeated indexes is implicitly adopted, and Equation (7) has an important physical interpretation.Define Then, Equation ( 7) is rewritten as Equation ( 9) has the same form as the continuity equation of fluid mechanics.In this way, C j is the j-th component of the vector of the probability current, so that Equation ( 9) expresses the conservation of the probability.From that it is deducted that the C j 's must be zero at the boundaries of the existence domain, which in most cases is infinite; hence, at the boundaries the JPDF p is zero.So, it is demonstrated that the FPK equation describes a convection problem.Now, for the clarity's sake some solutions of the reduced (steady-state) FPK equation are given.The motion equation of a second order oscillator with nonlinear restoring force is where W(t) is a unit strength Gaussian white noise.It is assumed that F(x) is a conservative force, which is it is the derivative of a potential function U(x) as Equation ( 10) is equivalent to the following two first order Itô's stochastic differential equations in the phase space: where . The FPK equation governing the stationary JPDF where 1 2 x x , and the current values of the variables have been denoted with lower case letters.The solution where C is a normalization constant.
In [5], the FPK equation is solved too for an N-degreeof-freedom (2N states) system, whose motion equation have the form where no summation has to be performed with respect to the index i; the white noises W i (t) are uncorrelated.The JPDF of the variables , In Equation ( 5) the equipartition of the energy among the degree of freedoms is notable.
Liu [8] considers a class of N-degree-of-freedom nonlinear discrete dynamic systems, whose equations of motions are where are constant parameters, and the Gaussian white noise processes W i are characterized by He gives the following solutions.I-Under the hypothesis:  the C j 's i are the same for all the degree of freedom;  the matrix ij b     is not singular so that the inverse matrix Δexists; the following relationships hold 2 where if , , and if , y y the probability flux is null around the existence domain.In this case.the steady state PDF is given by where the summation with respect to the indexes j and k is implicit.Equation ( 20) may or may not result in a closed form depending on whether the integral is analitically evaluable.
II-The addenda in Equation ( 9) are singularly zero according to the relationships where L j are first order differential operators, while g j and h j are nonlinear functions.Then, the JPDF is In [13] by applying the principle of detailed balance, the steady state JPDF's of some notable oscillators are found.The first of them is where   h  is a generic function of the mechanical energy of the oscillator, that is The steady-state JPDF is Consider the oscillator with parametric excitation where   h  is a function of the mechanical energy , and the processes W 1 , W 2 , W 3 are uncorrelated Gaussian white noises with power spectral densities 1 2 , , K K and 3 K , respectively.The steady-state JPDF is where   , x x


 is given by the equation A more general case of FPK equation for an oscillator with parametric excitation of white noises is solved in [15] by applying the concept of generalized stationary potential.It is where the summation rule with respect to a repeated index is valid, are nonlinear functions, and the white noises W i are characterized by the correlation . The steady-state PDF has the form of Equation ( 27) with where  is the generalized stationary potential, 2 1 2 y x   , and   a function of the mechanic energy.Equation ( 30) is valid under the condition    are the derivatives of  with respect to x, y, and to y twice, respectively; 0  is the derivative of  with respect to , that is the generic value of .
At this stage the functions D(x) and  0 () are still unknown: they can be identified by applying Equation (31); see the examples in [15].It is notable that Equation ( 31) is a very restrictive relation between the system nonlinearity and the excitation parameters, which is rarely met in practice.
Zhu found the exact stationary solution of the FPK equation pertaining to several classes of nonlinear dynamic systems [18,20,26,27].First, we recall the case of the following nonlinear stochastic system: In Equation ( 32) a and b are constants, 2 where the K ij 's define the correlations among the white noises, that is Another class of dynamic systems for which the associated FPK equation has been extensively studed is that of stochastically excited and dissipated Hamiltonian systems with n-degree-of-freedom [17,19,20,26,27]: where Q i and P i are generalized displacements and momenta, respetively; is a twice differentiable Hamiltonian; are differentiable functions; In the second of Equations ( 35) the excitations are multiplicative: hence, for transforming them into Itô type stochastic differential equations, they are to be modified by means of the Wong-Zakai-Stratonovich corrective terms [55,56].The transformed equations in incremental form are where in general Ĥ H  and .

ij ij
The FPK equation associated with Equation ( 36) is where The steady-state PDF solution to Equation (37) has the form where the vectors q and p contain the generalized displacements and momenta, respectively.If the conservative Hamiltonian system , ,  depends on these, and it is the solution of n first order PDE's.
For both cases the reader is referred to [20] for the details.

General Formulation
In the finite element (FE) procedure the state space is discretized into a number of finite regions: the finite element mesh consists of a grid of points at which the JPDF p(z,t) is to be computed.It must be emphasized that in general the JPDF p exists in an infinite hyperdomain.If the FE procedure is applied to an infinite domain, three different approaches are possible: 1) the inner region of the domain is divided into finite elements, while the outer region is discarded since the quantity to be computed is assumed negligible therein; 2) the inner region of the domain is divided into finite elements, while the outer region is modelled by infinite elements that extend to infinity; 3) the outer region is modelled using boundary elements or series solution.The first approach is used herein as the probability density functions decay to small or very small values, when the variables are larger than a few standard deviations.An estimate of the standard deviations of the response variables is obtained easily by applying the equivalent linearization method [61].Inside an element s (Figure 1) the JPDF p(z,t) is approximated by p (els) , which is given by In Equation ( 39) N no are the nodes of the element, where the JPDF assumes the nodal values   s k p .Clearly, inside an element the JPDF varies according to the shape functions N k (z).The nodal values of contiguous elements must be equal.
If linear shape functions are selected, the values of the JPDF at the nodes are the only as unknowns.On the other hand, linear shape functions yield C 0 continuity, while C 1 continuity is required as the FPK equation is of second order.This problem is overcome by considering a weak solution of it.Define the operators For any weighting function (z) it must be in which the dependence on z has been omitted, and  denotes the domain of existence of p(z,t).Integrating the right-hand-side of Equation ( 41) by parts, and accounting for the fact that as 0 p    z , we obtain the weak form of the FPK equation as Equation ( 42) is obtainable also with a variational formulation [48].
If Equation ( 39) is inserted in (42), and the integration over the domain  is replaced by the sum of the integrals over each element, we obtain the matrix system in which p(t) is an N no   column vector containing the nodal values, and 0 is an N no   column vector whose components are zero.The various entries of the N no  N no matrices C, K depend on the choice that has been made for the weighting functions (z).Bergman and Heinrich [45], Langtangen [48], Köylüoglu and collaborators [50] use weighting functions different from the shape functions, while the other authors use weighting functions equal to the shape functions.This choice, which is referred as Bubnov-Galerkin FE method, gives   , ( 45) where is a generic element of the local "stiffness" matrix, m is an element of the local matrix multiplying , and E denotes that the integration is performed over the element.Transformation of each element matrix into global coordinates and summation over the elements allows the construction of the matrices C and K.
If the shape and weighting functions are chosen to be different, Equations (44,45) are replaced by, respectively Open Access ENG C. FLORIS 981

Refinements of the Method
Computer codes have been written and implemented to solve both the transient and the reduced FPK equations (in this case Equation ( 43) reduces to Kp = 0): because of brevity's sake the features of these codes are not expounded.It is only recalled that Equation ( 43) is solved by adopting a forward finite difference scheme.When the shape and weighting functions are equal, in the case of a rectangular element (Figure 1) the functions in local coordinates ,   are This choice has the undoubted advantage of making the solutions of the integrals faster (obviously, the integrals are numerically evaluated).
in which the constants a, b depend on the transformation from global to local coordinates, and are assumed differently in the different cases.Thus, the i  's are multiplied by Gaussian like functions that cause them to decay fast from the nodes.
Computing the JPDF of the state variables is not the final stage of the analysis of a random dynamic system.In fact, one needs to obtain the response statistics such as marginal densities, statistical moments, and mean upcrossing rate (MUR) functions.The former can be directly obtained as geometrical sections of the hypersurface having p(z,t) as equation.However, the other quantities require integration: of the marginal PDF as regards the response moments, and of the JPDF of X, X  as regards the MUR function.
In many cases, accepting or rejecting a JPDF or a PDF obtained by means of the FE method is decided by a visual inspection.On the contrary, a sound quantitative criterion of acceptance must be based on the errors that the computed statistics have with respect to the exact or simulative ones.In this paper, an improvement of the method is proposed in order to compute better estimates of the statistics without increasing the number of nodes.The nodal values are interpolated by both quadratic and cubic splines, the former choice bringing to the well known Cavalieri-Simpson's rule of integration.

Applications
The present section is devoted to the presentation of the results that are obtained applying the FE method to the solution of the FPK equations of different stochastic differential equations.In order to limit the computing time, the transient state is analyzed only for two scalar systems.In fact, the solution of the complete Equation ( 43) is much more time consuming as it requires a multi-fold discretization in the state coordinate and in time.On the contrary, the solution of the steady-state equation Kp = 0 can be easily performed even on a normal personal computer.Two-state systems are analyzed in this case.The results of the Bubnov-Galerkin version of the FE method are compared with those of the Petrov-Galerkin one, which is indicated as WRM.In computing the staistical moments, the effets of the spline interpolation of the nodal values are highlited.

Scalar Systems
The feasibility of the FE computer code that has been implemented is firstly tested by determining the time evolution of the PDF of the solution X of two scalar systems excited by a stationary Gaussian white noise.It has been chosen the Langevin equation with linear or nonlinear drift term.In the first case the PDF of X is Gaussian.The equation of motion are in the two cases, respectively.The steady-state PDF of X in Equation ( 51) is in which C is a normalization constant.
The analyses have been accomplished with the following values of the parameters: in Equation (50); 51) in such a way that the PDF (52) is bimodal.In both cases the initial conditions are random, that is is a Gaussian PDF with zero mean and standard deviation .Since the two systems are scalar, only 32 finite elements suffice to get a good agreement among numerical and theoretical results, but the computations  50) and ( 51), respectively.It is not surprising that the PDF of X depicted in Figure 2 converges very fast to the steady-state one as this is Gaussian like the initial one.However, the convergence is fast as in the case of the nonlinear Langevin equation even if the initial and the final PDF are quite different.

Duffing Oscillator
The second case that has been examined regards a Duffing oscillator excited by a Gaussian white noise.The governing equation of motion is where  and  are constant parameters, and W(t) is a unit strength Gaussian white noise.In this case, the steadystate JPDF is known to be [5]   in which , and C is a normalization constant, whose expression is where  K  is the modified Bessel function of order 1/4.It is recalled that, when both  and  are positive, the PDF of X is unimodal and Gaussian like, while for    the PDF is bimodal.An analytical solution is known also for the MUR function: The computations have been performed for

Oscillator with Parametric Excitation
Now, consider the dynamic system , ( 57) in which W 1 (t) and W 2 (t) are Gaussian uncorrelated white noises with power spectral densities K 1 and K 2 , respec

Open Access ENG
C. FLORIS 984 tively.The oscillator of Equation (57), which has the parametric excitation W 1 (t), belongs to the class of generalized stationary potential [14] if the condition   4 0 1 2

K K
 is fulfilled.The steady-state JPDF has the form (27) with It is notable that  , x x   does not depend on the in-  tensity of the parametric excitation W 1 (t).The computations have been performed for the following set of parameters: In the Bubnov-Galerkin approach, meshes with 400 and 900 elements have been tested.In the Petrov-Galerkin approach, the elements over one quarter of the integration domain are 400, and the constants a, b of Equation ( 49) are both 0.150.
The exact stationary PDF of X and those obtained by using the FE method are compared in the top plot of Figure 6.Again, both methods and all levels of discretization give good results so that they cannot be discriminated by a visual inspection.It is not possible to distinguish the PDF of the Bubnov-Galerkin approach from that obtained with different weighting and shape functions.The bottom plot compares the MUR functions   X x


. In this case, an analytical expression does not exist: the results labeled as exact are obtained inserting the exact JPDF in the Rice's formula, then the integral is numerical evaluated.The same procedure is used to obtain


from the FE results.The agreement is quite satisfactory.
The mean square values of X and X  are in Table 3.There are reported only the values obtained by means of 900 elements in the Bubnov's approach and 400 in the Petrov's one.It is evident that the errors are small in any case.However, Petrov's procedure allows a non negligible computing time saving.

Oscillator with Cubic Damping and Stiffness
Consider the following oscillator     Because of space limitations only the marginal PDF's of X for  = 1 are shown in Figure 7, and only for 800 elements over one quarter of the domain with shape functions different from weighting functions.A perfect agreement between the two curves can be noted.This is confirmed by the statistical moments: by simulation, and by means of the FE method.However, the former requires more than 9 hours of elaboration on a workstation with quadcore processor, while the latter few minutes.

Conclusions
This paper aims at giving a contribution to some questions regarding the FE method for solving the FPK equatio n, which remains unanswered or partially answered.Firstly, the existence domain of a JPDF p(z,t) is generally infinite, while the FE solution is restricted to a limited domain.This fact has two shortcomings: 1) the former is the necessity of having a good estimate of the response standard deviations in order to give the integration do- main appropriate dimensions; 2) no precise information is obtained on the values of the JPDF in its tails, which are very important in structural reliability.An enlargement of the domain of integration causes a notable increase of the computations in many cases.According to Langley [47], a combination of finite elements with boundary elements or infinite elements does not seem to be effective in the case of the FPK equation.
A second question is in some way related to the former: spurious oscillations of the PDF values and meaningless regions of negative values are possible, particularly in the tails, when the mesh is too coarse or badly defined.Finally, the application of the FE method to problems with more than 3 dimensions is an open question because of computer storage and speed requirements.
Computer programs have been written and implemented for both formulations of the FE method: the Bubnov-Galerkin method, which uses equal weighting and shape functions, and the Petrov-Galerkin method in which these functions are different.
Examining the results that have obtained for different stochastic systems with additive or both additive and multiplicative excitations, the following conclusions can be drawn: 1) The FE method is very feasible for solving the reduced FPK equation when the response states are few, say not more than three or four.If there are more states, programming is more cumbersome, and the computing times are unavoidably longer, not competitive with Monte Carlo simulation.2) The steady-state statistical moments of the response are computed with small errors and a reasonable charge for computations.For two states, this is lesser than that required by a Monte Carlo simulation, and larger than that required by a moment equation approach.The FE method has the advantage that no closure schemes are required for computing the moments differently from the moment equation approach.
3) If one wants to limit the charge for computations, the Open Access ENG C. FLORIS 986 integration domain must be kept into 4 -5 standard deviations of the variables.In this way, it is not possible to get reliable values in the tails of the JPDF, in which negative values can be encountered.4) The use of weighting functions different from shape functions has given encouraging results so that further study in this direction seems to be promising.5) The FE method is proved to be feasible even for transient responses, but in this case the computing time increases dramatically as the equations are to be solved in every time instant.
0.01 s, that is t  0.01.The time evolutions of   , X p x t are in Figures 2 and 3 for Equations (

Figure 2 .
Figure 2. Time evolution of the PDF p(x) of X in Equation (50), random initial condition: stars FE solution, line steady-state Gaussian PDF.

Figure 3 .
Figure 3.Time evolution of the PDF p(x) of X in Equation (51), random initial condition: stars FE solution, line steady-state exact PDF.

Figure 4 .
Figure 4. Steady-state statistics of X for Duffing oscillator with  = 1: top rendering of the two-dimensional JPDF of X X  , 

Figure 5 .
Figure 5. Steady-state statistics of X for Duffing oscillator with  = 1: top rendering of the two-dimensional JPDF of X X  , ; middle marginal PDF of X; bottom MUR function   X x 

Figure 6 .
Figure 6.Steady-state statistics of X in Equation (57): top marginal PDF of X; bottom MUR function   X x 

1 
  , and ,  are nonlinear parameters.The FPK equation associated with this oscillator does not admit an analytical solution.Thus, the finite element estimates have been compared with those obtained by means of Monte Carlo simulation.Sets with 60,000 samples of motion histories have been simulated.The parameters take the following values: 0 2π,

Figure 7 .
Figure 7. Steady-state PDF of X in Equation (59): blue line WRM FE estimate, black line monte carlo simulation).
[60]ts first integral is the Hamiltonian itself[60], the function  depends on H only, and it is by solving a system of first-order linear ordinary differential equations.If there are n independent integrals of motion 1 obtained