Application of multi-scale finite element methods to ... - Semantic Scholar
Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526 www.elsevier.com/locate/cma
Application of multi-scale finite element methods to the solution of the Fokker–Planck equation Arif Masud a
a,1
, Lawrence A. Bergman
b,*
Department of Civil and Materials Engineering, University of Illinois at Chicago, M/C 246, 842 W. Taylor St. Chicago, IL 60607, USA b Department of Aeronautical and Astronautical Engineering, University of Illinois at Urbana-Champaign, 104 S. Wright St., 321E Talbot Lab., MC-236, Urbana, IL 61801, USA Received 21 November 2003; received in revised form 10 May 2004; accepted 9 June 2004
Abstract This paper presents an application of multi-scale finite element methods to the solution of the multi-dimensional Fokker–Planck equation. The Fokker–Planck, or forward Kolmogorov, equation is a degenerate convective–diffusion equation arising in Markov-Process theory. It governs the evolution of the transition probability density function of the response of a broad class of dynamical systems driven by Gaussian noise, and completely describes the response process. Analytical solutions for the Fokker–Planck equation have been developed for only a limited number of lowdimensional systems, leading to a large body of approximation theory. One such approach successfully applied to the solution of these problems in the past is the finite element method, though for systems of dimension three or less. In this paper, a multi-scale finite element method is applied to the Fokker–Planck equation in an effort to develop a formulation that can yield higher accuracy on cruder spatial discretizations, thus reducing the computational overhead associated with large scale problems that arise in higher dimensions. 2004 Elsevier B.V. All rights reserved. Keywords: Multi-scale finite element methods; Fokker–Planck equation
Background Characterization of the response of dynamical systems driven by random noise has been a subject of continuing interest for more than 60 years. An extensive body of literature on the subject exists, and the interested reader can consult one or more recent review articles (e.g., [1]) for an overview. The focus here, *
Corresponding author. Tel.: +1 217 333 4970; fax: +1 217 244 0720. E-mail addresses: [email protected] (A. Masud), [email protected] (L.A. Bergman). 1 Tel.: +1 312 996 4887; fax: +1 312 996 2426. 0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.06.041
1514
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
though, will be on one particular aspect of the problem: numerical solution of the Fokker–Planck equation. It is well known that the response of a dynamical system, without memory, driven by Gaussian white noise, forms a Markov process. The evolution of the response, then, can be posed as an integral equation, the Chapman–Kolmogorov–Smulochowski equation, or its partial differential equation equivalents, the forward (Fokker–Planck) and backward Kolmogorov equations. For a dynamical system of dimension 2n, the Fokker–Planck equation is an initial-boundary value problem in 2n spatial dimensions plus the time dimension. The equation is convective–diffusive and non-self-adjoint. In general, analytical solutions exist, at least in terms of quadratures, only for scalar dynamical systems. There has been some success with systems of dimension two, particularly in determining stationary, or long term, behavior. This absence of analytical results has motivated the development of a body of approximation theory for estimating the probability density function of response, which includes but is not restricted to maximum entropy distributions, cell mapping, path integration, Monte Carlo methods, and numerical solution of the Fokker–Planck equation, all of which are discussed in some detail, with references, in Schueller [1]. This paper focuses on the last of these and presents an advantageous new finite element procedure for the solution of the Fokker– Planck equation that can be applied to problems of higher dimension. The earliest work on numerical solution of the Fokker–Planck equation applied to realistic engineering systems of dimension two and three appears to have been published by Atkinson [2] and Wen [3,4]. Both solved Fokker–Planck equations associated with non-linear systems using a global Galerkin-type approach. Bergman and Heinrich [5] reported the first successful application of the finite element method to the solution of the related backward-Kolmogorov equation, the formal adjoint of the Fokker–Planck equation, to solve the first passage problem for the second order linear system driven by additive Gaussian white noise; Bergman and Spencer [6,7] applied the finite element method to the solution of the Fokker–Planck equation for a number of two dimensional systems subjected to additive white noise; and Wojtkiewicz et al. [8] examined several three dimensional systems. More recently, Wedig and Wagner [9] examined stationary solutions of systems of dimension four and six, using a global Galerkin-based formulation. Our objective in this work is to develop a new method that possesses higher accuracy on crude discretizations. This would substantially lower the computational costs associated with larger meshes that arise in the analysis of large scale dynamical systems that are subjected to additive and multiplicative Gaussian excitations, retaining a sufficient number of states to capture all of the important behavior. To that end, we present the development of a multi-scale finite element method and its application to the solution of the Fokker–Planck equation and examine its efficacy through model applications.
1. Introduction We seek to develop a stabilized formulation in higher dimensions so that we can apply numerical solution techniques such as the finite element method to the Fokker–Planck equation or to its formal adjoint, the backward Kolmogorov equation. The Fokker–Planck, or forward Kolmogorov, equation is a degenerate convective–diffusion equation arising in Markov-process theory. The first applications of the stabilized methods to the convective–diffusion equation can be traced back to the early 80s when Hughes and colleagues [10,11] introduced the Streamline-Upwind-Petrov–Galerkin (SUPG) method to address the issue of lack of stability of the standard Galerkin approach to this problem. The SUPG method turned out to be the forerunner of a new class of stabilization schemes, namely the Galerkin/Least-Squares (GLS) stabilization methods [12]. In the GLS method a least-squares form of the residuals, based on the corresponding Euler–Lagrange equations, is added to the Galerkin finite element formulation. The underlying philosophy of the stabilized methods is to strengthen the classical variational formulations so that discrete approximations, which would otherwise be unstable, become stable and convergent. In the mid 90s Hughes [13] revisited the origins of the stabilization schemes from a variational view point and presented the variational
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1515
multi-scale method. In this method the different stabilization techniques come together as special cases of the underlying sub-grid scale modeling concept. Using this method Masud and Hughes [14] developed a new stabilized formulation for the Darcy equation and showed the superior performance of the method as compared to the GLS method. The multi-scale decomposition of the underlying fields was employed by Masud and coworkers to develop stabilized formulations for the incompressible Navier–Stokes equations [15], the convective–diffusive heat transfer problem [16], and the advection–diffusion equation [17]. We approach the present topic with a new point of view that finds its roots in Hughes Variational Multiscale method. We know that all numerical methods are based on the notion of the existence of underlying computational grids (i.e., a discrete set of points on which the field variables are approximated). The distance between these points is the characteristic length scale of discretization and is represented by h. This length scale, which is strictly an artifact of discretization, plays a crucial role in the computed solution. For a well-defined problem and with a good numerical procedure, one is able to resolve the features of the solution that are larger than, or at least of the order of, the discretization scale h. However, features of the solution which are finer than h are not captured. Within the context of Galerkin finite element methods, the only way to capture these fine scale features is to refine the mesh; i.e., to reduce the artificial length scale h. It is important to realize that if this approach is adopted, one very quickly reaches the limit of any current computational platform. Thus, we consider here whether it is possible to capture all features of the solution, both coarse and fine, on coarse meshes. In other words we endeavor to develop a numerical scheme that shows high accuracy on crude discretizations in an effort to reduce the computational overhead associated with large scale problems that arise in higher dimensions. With this in view, let us assume that the total solution f is decomposable into f ¼ f h þ f e; where f h is the part of the solution that can be obtained with a good numerical scheme and f e is the part of the solution that is lost because it is finer than the artificial length scale h. In other words, f e represents the error in the solution. We now present a systematic procedure wherein we build this part back into the formulation. The beauty of the procedure is that it automatically yields a stabilized form of the problem. In addition, we get an explicit expression for the stabilization parameter. It is important to note that since a stabilization term is not being a priori designed for the problem, the method is amenable to approximation in n-dimensional space. Thus, it is ideally suited to Kolmogorov equations in higher dimensions, where generalization of existing stabilization methods is not straightforward and often not effective. 2. The Fokker–Planck equation The Fokker–Planck, or forward Kolmogorov, equation corresponding to an n-dimensional dynamical system subjected to Gaussian white noise excitation is given by ! n n X m X of o 1 X o ¼ ðzj f Þ þ ½H ij f ; ð1Þ ot oxj 2 i¼1 j¼1 oxi oxj j¼1 where f(x, tjx0, t0) is the transition probability density function, vector x represents an n dimensional space, and z(x) 2 Rn, H(x) 2 Q Rn · Rm are the drift vector and diffusion matrix, respectively. The initial condition is n given by f ðx; 0jx0 Þ ¼ i¼1 dðxi x0i Þ. Let def
r ¼ gradient operator in ‘n’ dimensions; def
r ¼ divergence operator in ‘n’ dimensions; def
D ¼ Laplace operator in ‘n’ dimensions:
1516
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
Furthermore, let X Rn be an open bounded region in n dimensions with piecewise smooth boundary C. We assume n P 2. The Fokker–Planck equation can be written in direct notation as f_ ¼ r ðzf Þ þ 12DðHf Þ in X 0; T ½:
ð2Þ
Eq. (2) together with the initial condition constitutes the strong form of the problem. Now, let w H1(X) be the weighting function for the probability density function. Taking the inner product of the weighting function with (2) and applying the divergence theorem yields the weak form of the problem. ðw; f_ Þ þ ðw; r ðzf ÞÞ þ 12ðrw; Hrf Þ ¼ 0 in X 0; T ½;
ð3Þ
where (Æ, Æ) is the L2 inner product. In the final step, we let V H 1 ðXÞ \ C 0 ðXÞ represent the space of piecewise continuous trial solutions and weighting functions. Substituting the discrete counterparts of f and w into (3) yields the standard Galerkin form of the problem. Remark. In general, the standard Galerkin form requires a stabilized method to yield a solution. We first present the salient features of the multi-scale decomposition that we will employ in the present development. The Variational Multiscale Method was first proposed by Hughes [13].
3. Multiscale decomposition We consider the bounded domain X discretized into numel non-overlapping regions Xe (element domains in n dimensions) with boundaries Ce, e = 1, 2, . . . , numel such that X¼
n[ umel
Xe :
e¼1
We denote the union of element interiors and element boundaries by X 0 and C 0 , respectively, X0 ¼
n[ umel
ðintÞXe
ðelement interiorsÞ;
ðintÞCe
ðelement boundariesÞ
e¼1
C0 ¼
n[ umel e¼1
and we assume an overlapping sum decomposition of the PD (probability density) field into coarse, or resolvable-scales and fine or unresolved-scales. f ðxÞ ¼ f ðxÞ þ f 0 ðxÞ : ð4Þ |ffl{zffl} |ffl{zffl} coarse scale
fine scale
Likewise, we assume an overlapping sum decomposition of the weighting function into coarse and fine and w 0 , respectively scale components, indicated by w ðxÞ þ w0 ðxÞ : wðxÞ ¼ w |ffl{zffl} |fflffl{zfflffl} coarse scale
ð5Þ
fine scale
We further assume that the fine scales, although non-zero within the elements, vanish identically over the element boundaries
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
f 0 ¼ w0 ¼ 0
on C0 :
1517
ð6Þ
We specify a direct sum decomposition on the spaces of functions for the coarse and fine scale fields V0 ; V¼V
ð7Þ
in (7) is the space of trial solutions and weighting functions for the coarse scale PD field, and is where V identified with the standard finite element space C 0 ðXÞ \ V; f 2 V
e Þ ¼ Pk ðXe Þ: VðX
Here, Pk ðXe Þ denotes the complete polynomials of order k over Xe. On the other hand, various characterizations of V0 are possible, subject to the restriction imposed by the and V0 to be linearly independent. Consequently, in the discrete stability of the formulation that requires V 0 case, V can contain various finite dimensional approximations; e.g., bubble functions or p-refinements that satisfy Eq. (6), i.e., f 0 2 V 0 ¼ ff 0 jf 0 ¼ 0 on C0 g:
ð8Þ
4. The multiscale problem We now substitute the trial solution (4) and the weighting function (5) in the standard weak form (3). This becomes the point of departure from the conventional Galerkin formulation 0 ð w þ w0 ; f_ þ f_ Þ þ ð w þ w0 ; r ðzðf þ f 0 ÞÞÞ þ 12ðrð w þ w0 Þ; Hrðf þ f 0 ÞÞ ¼ 0:
ð9Þ
Remarks (1) With suitable assumptions on the fine scale field, as stipulated in (6), and employing the linearity of the weighting function slot, we can split the problem into coarse and fine scale sub-problems, indicated as and W0 , respectively. W (2) We assume that f 0 is represented by piecewise polynomials of sufficiently high order, continuous in x but discontinuous in time. In particular, f 0 is assumed to be composed of piecewise constant-in-time functions [15]. Therefore f_ ¼ f_ . The coarse scale sub-problem, W ð w; f_ Þ þ ð w; r ðzðf þ f 0 ÞÞÞ þ 12ðr w; Hrðf þ f 0 ÞÞ ¼ 0:
ð10Þ
contains only the coarse scale components. The weighting function slot in W The fine scale sub-problem, W0 ð w; f_ Þ þ ð w; r ðzðf þ f 0 ÞÞÞ þ 12ðrw0 ; Hrðf þ f 0 ÞÞ ¼ 0:
ð11Þ
The weighting function slot in W0 contains only the fine scale components. The general idea at this point is to solve the fine scale problem to obtain the fine scale solution f 0 which can then be substituted into (10), thereby eliminating the fine scales yet retaining their effect.
1518
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
5. Solution to the multi-scale problem 5.1. Solution to the fine scale problem Consider now the fine scale part of the weak form W0 , which because of the assumption on the fine scale space, is defined over the sum of element interiors indicated as X 0 . Exploiting linearity of the solution slot on the left hand side of (11), we have 1 1 ðw0 ; f_ ÞX0 þ ðw0 ; r ðzf ÞÞX0 þ ðw0 ; r ðzf 0 ÞÞX0 þ ðrw0 ; Hrf ÞX0 þ ðrw0 ; Hrf 0 ÞX0 ¼ 0: |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 2 2
ð12Þ
Term 3
Further, consider Term 3 alone in Eq. (12) and apply the mathematical identity r ðzf 0 Þ ¼ f 0 div z þ z rf 0 ; so that ðw0 ; r ðzf 0 ÞÞX0 ¼ ðw0 ; f 0 div zÞX0 þ ðw0 ; z rf 0 ÞX0 :
ð13Þ
Now, apply integration by parts to the left-hand side of (13) ðw0 ; ðz nÞf 0 ÞC0 ðz rw0 ; f 0 Þ ¼ ðw0 ; f 0 div zÞX0 þ ðw0 ; z rf 0 ÞX0 : 0
ð14Þ 0
In (14), the first term on the left-hand side is zero because w = 0 on C (see Eq. (6)). For a piece-wise constant projection of z div z ¼ 0: Accordingly, the first term on the right-hand side is zero, and we see that the right-hand side is equal to the negative of the left-hand side. (Note: Switch positions of w 0 and f 0 to see this.) This proves that the third term in Eq. (12) is zero. Accordingly, rearranging terms in (12) yields 1 ðrw0 ; Hrf 0 ÞX0 2
¼ ðw0 ; f_ ÞX0 ðw0 ; r ðzf ÞÞX0 12 ðrw0 ; Hrf ÞX0 ¼ ðw0 ; f_ þ r ðzf Þ 12DðH f Þ ÞX0 : |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
ð15Þ
¼r
Remark. Consider the bracketed term on the right-hand side in (15) and compare it with (2). It is clear that this is the residual of the coarse scales for the Euler–Lagrange equation of the problem i.e., the fine scales are in fact driven by the residual of the coarse scales. Having now shown this, our objective is to solve (15) to get an expression for f 0 that can be substituted into (10) (the coarse scale problem), thereby eliminating f 0 yet retaining its effect. Without loss of generality, we assume that the fine scales are represented via bubble functions over the element domain; i.e., f 0 jXe ¼ be be
on Xe ;
ð16Þ
w0 jXe ¼ be ce
on Xe ;
ð17Þ
where be represents the bubble shape function over element domains, and be and ce represent the coefficients for the fine scale trial solutions and weighting functions, respectively. Substituting (16) and (17) into the fine
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1519
scale problem (15), taking the vectors of constant coefficients be and ce out of the integral expressions, and employing the arbitrariness of ce results in Z 2 e e b ¼R b r dX : ð18Þ T ðrbe Hrbe Þ dX X0 X0 We can now reconstruct the fine scale field over the element domain using (16). " Z # 2 f 0 ðxÞ ¼ be R ber dX : eT e 0 X 0 ðrb Hrb Þ dX X
ð19Þ
For multi-linear elements, the spacial part of the residual is constant. To simplify the expression we consider the constant projection of the f_ term. Accordingly, the fine scale field f 0 (x) can be written in a concise form as f 0 ðxÞ ¼ sr; where s¼R
2be X0
R
X0 eT
ð20Þ be dX
ðrb Hrbe Þ dX
ð21Þ
:
It is important to note that the matrix H possesses only one non-zero entry, i.e., the term in the hmm slot. Therefore, the final form of the stabilization parameters for the Fokker–Planck equation is s¼R
R 2be X0 be dX : ðbe;m jhmm jbe;m Þ dX X0
ð22Þ
Remarks (1) It can be proven that stabilization parameter s is a positive bounded function. (2) It is important to note that the choice of the bubble function only affects the definition of the stabilization parameter s. Given the expression for f 0 (the fine scale field), we next substitute it into (10) and get the modified form for W. 5.2. Solution to the coarse scale problem Starting with (10) and employing the linearity of the solution slot, we expand to get ð w; f_ Þ þ ðð w; r ðzf ÞÞ þ ð w; r ðzf 0 ÞÞ þ 12ðr w; Hrf Þ þ 12ðr w; Hrf 0 ÞÞ ¼ 0:
ð23Þ
Applying integration by parts to the third and fifth terms, and applying (6) gives 1 1 Þ; f 0 Þ ¼ 0: w; Hrf Þ ðDðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ ðr ðz wÞ; f 0 Þ þ ðr |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 2 2 |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} Term 3
Term 5
ð24Þ
1520
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
Combining terms 3 and 5 and substituting f 0 from (20) in (24) Þ; srÞ ¼ 0: w; Hrf Þ þ ðr ðz wÞ þ 12DðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ þ 12ðr Substituting r from (16) yields the modified formulation Þ; s f_ þ r ðzf Þ 12DðH f Þ ¼ 0: w; Hrf Þ þ r ðz wÞ þ 12DðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ þ 12ðr
ð25Þ
ð26Þ
Remarks (1) The first three terms in (26) are the Galerkin terms. (2) The fourth term in (26) has appeared because of the assumption of the existence of fine scales in the problem. This term is not present in the classical Galerkin formulation for the Fokker–Planck equation, and represents f 0 in the problem. (3) The fourth term is in fact a residual-based term; i.e., the trial solution slot contains the residual of the Euler–Lagrange equations for the coarse scales. Accordingly, the formulation is consistent and accommodates the exact solution. (4) The additional (fourth) term enhances the stability of the formulation, without upsetting the consistency. It is important to note that stability and consistency lead to convergence. (5) The consistent mass matrix in this formulation gets an additional contribution from the stabilization term Þ; sf_ : ð w; f_ Þ þ r ðz wÞ þ 12DðH w ð27Þ (6) For the case of linear elements, a comparison of the modified formulation with the two celebrated staÞ bilized methods, i.e., the SUPG and the GLS methods is straightforward. If we drop the þ 12 DðH w term from the weighting function slot of the stabilization term, we recover the SUPG formulation Þ term in the weighting function for the Fokker–Planck equation. If we change the sign of þ 12 DðH w slot of the additional stabilization term and add the weighting counterpart for the rate term, we recover the GLS formulation for the Fokker–Planck equation. However, it is important to note that the stabilization terms in SUPG and GLS contain stabilization parameters that designed based on the consideration of achieving nodally exact solution in the one-dimensional setting. Generalization of the nodally exact stabilization parameters in 1-D to the case of nodally exact parameters in multi-dimensions has been a challenging task. A literature review in the field of fluid mechanics would reveal various attempts in designing the optimal stabilization parameters. (7) In the proposed formulation the definition of the stabilization parameter s appears naturally and is given in Eq. (21). This is in sharp contrast to the SUPG or the GLS stabilization schemes for the Kolmogorov equations. (8) It is important to note that the use of bubble functions resides completely in the definition of s, as given in Eq. (22). (9) The modified formulation is completely expressed in terms of the coarse scale PD function f . This means that the dimension of the element stiffness matrix is the same as that for the Galerkin method. (10) The proposed formulation is stable in Rn, where n is the spatial dimension of the Fokker–Planck equation considered. (11) The formulation accommodates Lagrange-based C0 functions. In compact notation, the multi-linear shape function can be expressed as
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
N a ðnÞ ¼
n 1 Y ð1 þ nai ni Þ 2n i¼1
1521
ð28Þ
n ¼ n1 ; n2 ; . . . nn : (12) The formulation also accommodates triangles, tetrahedron and their generalizations in n-dimensions.
6. Numerical examples We have developed 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals (see Fig. 1) for application to the Fokker–Planck equation. We have employed the simplest polynomial bubbles, which in the element natural coordinate frame are as follows: be ðn; gÞquadrilateral ¼ ð1 n2 Þð1 g2 Þ;
ð29Þ
be ðn; gÞtriangle ¼ 27ngð1 n gÞ:
ð30Þ
6.1. The rate of convergence study The first numerical simulation is a study of convergence rates on a model problem. The domain under consideration is a bi-unit square with origin at (x, y) = (0, 0). The exact solution is assumed to be of the form f ¼ sin
2px 2py sin : L L
ð31Þ
The ‘‘forcing function’’ as evaluated from Eq. (1) is given by 2p 2px 2py 2px 2py p 2px 2py z1 cos sin z2 sin cos sin sin f_ ¼ : L L L L L L L L
Fig. 1. A family of 2-D linear and quadratic elements.
ð32Þ
1522
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
In specifying the boundary-value problem, the forcing function is integrated over X while f = 0 is prescribed at the boundary. Note that H22 = 1 in this convergence rate study. The meshes employed in the convergence study are regular structured meshes. For linear quadrilateral elements, the meshes consist of 100, 400, 1600 and 6400 elements. The linear triangular element meshes consist of exactly twice as many elements and are constructed by bisecting the quadrilaterals along the main diagonal. For quadratic quadrilateral elements, the meshes employed consist of 100, 400, 1600 and 6400 elements. Again, the quadratic triangular element meshes consist of twice as many elements. The element mesh parameter, h, is taken to be the edge length of the elements for the quadrilaterals, and the short-edge length for triangles. Convergence rates for jzij = 10 are presented in Figs. 2 and 3. Fig. 2 shows the L2-norm and H1-seminorm of the scalar field for the linear quadrilateral and triangle. We see that optimal rates are achieved. Fig. 3 presents the corresponding rates for the quadratic elements. Once again we see optimal convergence rates. Fig. 4a and b present the elevation plots of the exact and the computed fields for a representative
Fig. 2. Convergence rates for linear elements.
Fig. 3. Convergence rates for quadratic elements.
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1523
Fig. 4. (a) Elevation plot of the scalar field (exact solution). (b) Elevation plot of the scalar field for the 4-node quadrilaterals (20 · 20 element mesh).
mesh of 4-node elements. These plots present a qualitative comparison of the computed solution with the corresponding exact solution, and show a stable response of the elements. 6.2. The linear oscillator (2D calculations) Simple verification of the proposed formulation (26) is afforded by the solution of the Fokker–Planck equation for the two-dimensional linear oscillator subjected to additive white noise excitation, for which an exact solution exists. The non-dimensionalized Fokker–Planck equation is given by 1 of o2 f o o ¼ 2n 2 ð2ny þ xÞf ðyf Þ: 2p os oy oy ox
ð33Þ
Here, x = x1/r, y = x2/rx0, s = x0t/2p, and r2 ¼ pS 0 =2nx30 . The damping parameter is n = 0.2. The computational domain is given by 10 6 x, y 6 10; the extent of the domain, 10, is chosen based on experience as
0.16 0.14
fx (0,0,τ)
0.12 0.10 0.08 0.06 0.04 0.02 Present Calculations 0.00 0
2
4
6
8
10
12
14
16
18
time Fig. 5. Probability density function for the linear oscillator at (0, 0).
1524
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
necessary in transient computations with non-zero initial conditions to minimize the incidence of probability mass reflecting off the boundaries and back into the interior. The mesh employed in the simulation is composed of 50 · 50 4-node elements, and the time-integration method is the Crank–Nicholson algorithm with Dt = 0.005. Boundary conditions are set such that f = 0 at the boundaries. The initial condition is a Gaussian distribution with standard deviation of unity, centered at (x, y) = (3, 3). Fig. 5 shows the evolution of the probability density function at (0, 0). Note that, at stationarity, the value of the distribution at the origin converges to f(0, 0) = 1/2p, the exact solution. Fig. 6 shows the elevation plot of the probability density function at times 1, 2 and 4 s, superimposed on common axes. Fig. 7 shows the exact probability density function at stationarity, while Fig. 8 shows the corresponding computational result. There is no observable difference between the two to the machine precision.
Fig. 6. Evolution of the probability density function for the linear oscillator.
Fig. 7. Exact probability density function for the linear oscillator.
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1525
Fig. 8. Computed converged probability density function for the linear oscillator.
7. Conclusions We have developed a multi-scale formulation for the Fokker–Planck equation of a class of oscillators in Rn. The key point of the method is a multi-scale decomposition of the field variable into resolvable/coarse and unresolvable/fine scales. The proposed multi-scale method leads to a stabilized formulation which is free of any user-defined parameters. Since the fine or the subgrid scales in the formulation are proportional to the grid scale residual, the method is consistent. We have tested the formulation on linear and quadratic triangular and quadrilateral elements and have obtained optimum numerical convergence rates in the norms considered. The formulation was applied to a model problem, the two-dimensional linear oscillator subjected to external white noise excitation, for which an exact solution exists. The computational result is a converged solution on a mesh of 2500 four node quadrilateral elements, representing an order of magnitude savings over previously reported converged results. While the two-dimensional result is not overly compelling, it does show significant promise for higher dimensional calculations. Prior FE solutions required on the order of 100 elements per spatial dimension for nominal convergence of a well-behaved dynamical system. Thus, a four-dimension problem which typically required 100 million elements would, incorporating the multiscale formulation, require nearly two orders of magnitude fewer. For the moment, efforts continue toward further refinement of the procedure with eventual application to larger scale systems subjected to both additive and multiplicative excitations.
References [1] [2] [3] [4] [5]
G.I. Schueller (Ed.), A state-of-the-art report on computational stochastic mechanics, Prob. Engrg. Mech. 13/4 (1997) 197–321. J.D. Atkinson, Eigenfunction expansions for randomly excited non-linear systems, J. Sound Vib. 30 (2) (1973) 153–172. Y.K. Wen, Approximate method for nonlinear random vibration, J. Engrg. Mech. Div., ASCE 101 (EM4) (1975) 389–401. Y.K. Wen, Method for random vibration of hysteretic systems, J. Engrg. Mech. Div., ASCE 102 (EM2) (1976) 249–263. L.A. Bergman, J.C. Heinrich, On the reliability of the linear oscillator and systems of coupled oscillators, Int. J. Numer. Methods Engrg. 18 (1982) 1271–1295. [6] L.A. Bergman, B.F. Spencer Jr., Robust numerical solution of the transient Fokker–Planck equation for nonlinear dynamical systems, in: N. Bellomo, F. Casciati (Eds.), Nonlinear Stochastic Mechanics, Springer-Verlag, Berlin, 1992, pp. 49–60. [7] B.F. Spencer, L.A. Bergman, Numerical solutions of the Fokker–Planck equations for nonlinear stochastic systems, Nonlinear Dyn. 4 (1993) 357–372.
1526
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
[8] S.F. Wojtkiewicz, L.A. Bergman, B.F. Spencer, Numerical solution of some three-state random vibration problems, in: S.C. Sinha (Ed.), Vibration of Nonlinear, Random, and Time-Varying Systems, ASME DE-84-1, 1995, pp. 939–947. [9] W.V. Wedig, Utz Von Wagner, On the calculation of stationary solutions of multi-dimensional Fokker–Planck equations by orthogonal functions, Nonlinear Dyn. 21 (2000) 289–306. [10] T.J.R. Hughes, A.N. Brooks, Streamline-upwind/Petrov–Galerkin methods for advection dominated flows, in: Proceedings of the Third International Conference on Finite Element Methods in Fluid Flow, Banff, June 1980, pp. 283–292. [11] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1–3) (1982) 199–259. [12] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (2) (1989) 173–189. [13] T.J.R. Hughes, Multiscale phenomena: Greens functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. [14] A. Masud, T.J.R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg. 191 (39–40) (2002) 4341–4370. [15] A. Masud, On a stabilized finite element formulation for incompressible Navier–Stokes equations, in: Proceedings of the Fourth US–Japan Conference on Computational Fluid Dynamics, Tokyo, Japan, May 2002. [16] M. Ayub, A. Masud, A new stabilized formulation for convective–diffusive heat transfer, Numer. Heat Transfer, Part B 44 (2003) 1–23. [17] A. Masud, R. Khurram, A multiscale/stabilized finite element method for the advection–diffusion equation, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018.
Application of multi-scale finite element methods to the solution of the Fokker–Planck equation Arif Masud a
a,1
, Lawrence A. Bergman
b,*
Department of Civil and Materials Engineering, University of Illinois at Chicago, M/C 246, 842 W. Taylor St. Chicago, IL 60607, USA b Department of Aeronautical and Astronautical Engineering, University of Illinois at Urbana-Champaign, 104 S. Wright St., 321E Talbot Lab., MC-236, Urbana, IL 61801, USA Received 21 November 2003; received in revised form 10 May 2004; accepted 9 June 2004
Abstract This paper presents an application of multi-scale finite element methods to the solution of the multi-dimensional Fokker–Planck equation. The Fokker–Planck, or forward Kolmogorov, equation is a degenerate convective–diffusion equation arising in Markov-Process theory. It governs the evolution of the transition probability density function of the response of a broad class of dynamical systems driven by Gaussian noise, and completely describes the response process. Analytical solutions for the Fokker–Planck equation have been developed for only a limited number of lowdimensional systems, leading to a large body of approximation theory. One such approach successfully applied to the solution of these problems in the past is the finite element method, though for systems of dimension three or less. In this paper, a multi-scale finite element method is applied to the Fokker–Planck equation in an effort to develop a formulation that can yield higher accuracy on cruder spatial discretizations, thus reducing the computational overhead associated with large scale problems that arise in higher dimensions. 2004 Elsevier B.V. All rights reserved. Keywords: Multi-scale finite element methods; Fokker–Planck equation
Background Characterization of the response of dynamical systems driven by random noise has been a subject of continuing interest for more than 60 years. An extensive body of literature on the subject exists, and the interested reader can consult one or more recent review articles (e.g., [1]) for an overview. The focus here, *
Corresponding author. Tel.: +1 217 333 4970; fax: +1 217 244 0720. E-mail addresses: [email protected] (A. Masud), [email protected] (L.A. Bergman). 1 Tel.: +1 312 996 4887; fax: +1 312 996 2426. 0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.06.041
1514
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
though, will be on one particular aspect of the problem: numerical solution of the Fokker–Planck equation. It is well known that the response of a dynamical system, without memory, driven by Gaussian white noise, forms a Markov process. The evolution of the response, then, can be posed as an integral equation, the Chapman–Kolmogorov–Smulochowski equation, or its partial differential equation equivalents, the forward (Fokker–Planck) and backward Kolmogorov equations. For a dynamical system of dimension 2n, the Fokker–Planck equation is an initial-boundary value problem in 2n spatial dimensions plus the time dimension. The equation is convective–diffusive and non-self-adjoint. In general, analytical solutions exist, at least in terms of quadratures, only for scalar dynamical systems. There has been some success with systems of dimension two, particularly in determining stationary, or long term, behavior. This absence of analytical results has motivated the development of a body of approximation theory for estimating the probability density function of response, which includes but is not restricted to maximum entropy distributions, cell mapping, path integration, Monte Carlo methods, and numerical solution of the Fokker–Planck equation, all of which are discussed in some detail, with references, in Schueller [1]. This paper focuses on the last of these and presents an advantageous new finite element procedure for the solution of the Fokker– Planck equation that can be applied to problems of higher dimension. The earliest work on numerical solution of the Fokker–Planck equation applied to realistic engineering systems of dimension two and three appears to have been published by Atkinson [2] and Wen [3,4]. Both solved Fokker–Planck equations associated with non-linear systems using a global Galerkin-type approach. Bergman and Heinrich [5] reported the first successful application of the finite element method to the solution of the related backward-Kolmogorov equation, the formal adjoint of the Fokker–Planck equation, to solve the first passage problem for the second order linear system driven by additive Gaussian white noise; Bergman and Spencer [6,7] applied the finite element method to the solution of the Fokker–Planck equation for a number of two dimensional systems subjected to additive white noise; and Wojtkiewicz et al. [8] examined several three dimensional systems. More recently, Wedig and Wagner [9] examined stationary solutions of systems of dimension four and six, using a global Galerkin-based formulation. Our objective in this work is to develop a new method that possesses higher accuracy on crude discretizations. This would substantially lower the computational costs associated with larger meshes that arise in the analysis of large scale dynamical systems that are subjected to additive and multiplicative Gaussian excitations, retaining a sufficient number of states to capture all of the important behavior. To that end, we present the development of a multi-scale finite element method and its application to the solution of the Fokker–Planck equation and examine its efficacy through model applications.
1. Introduction We seek to develop a stabilized formulation in higher dimensions so that we can apply numerical solution techniques such as the finite element method to the Fokker–Planck equation or to its formal adjoint, the backward Kolmogorov equation. The Fokker–Planck, or forward Kolmogorov, equation is a degenerate convective–diffusion equation arising in Markov-process theory. The first applications of the stabilized methods to the convective–diffusion equation can be traced back to the early 80s when Hughes and colleagues [10,11] introduced the Streamline-Upwind-Petrov–Galerkin (SUPG) method to address the issue of lack of stability of the standard Galerkin approach to this problem. The SUPG method turned out to be the forerunner of a new class of stabilization schemes, namely the Galerkin/Least-Squares (GLS) stabilization methods [12]. In the GLS method a least-squares form of the residuals, based on the corresponding Euler–Lagrange equations, is added to the Galerkin finite element formulation. The underlying philosophy of the stabilized methods is to strengthen the classical variational formulations so that discrete approximations, which would otherwise be unstable, become stable and convergent. In the mid 90s Hughes [13] revisited the origins of the stabilization schemes from a variational view point and presented the variational
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1515
multi-scale method. In this method the different stabilization techniques come together as special cases of the underlying sub-grid scale modeling concept. Using this method Masud and Hughes [14] developed a new stabilized formulation for the Darcy equation and showed the superior performance of the method as compared to the GLS method. The multi-scale decomposition of the underlying fields was employed by Masud and coworkers to develop stabilized formulations for the incompressible Navier–Stokes equations [15], the convective–diffusive heat transfer problem [16], and the advection–diffusion equation [17]. We approach the present topic with a new point of view that finds its roots in Hughes Variational Multiscale method. We know that all numerical methods are based on the notion of the existence of underlying computational grids (i.e., a discrete set of points on which the field variables are approximated). The distance between these points is the characteristic length scale of discretization and is represented by h. This length scale, which is strictly an artifact of discretization, plays a crucial role in the computed solution. For a well-defined problem and with a good numerical procedure, one is able to resolve the features of the solution that are larger than, or at least of the order of, the discretization scale h. However, features of the solution which are finer than h are not captured. Within the context of Galerkin finite element methods, the only way to capture these fine scale features is to refine the mesh; i.e., to reduce the artificial length scale h. It is important to realize that if this approach is adopted, one very quickly reaches the limit of any current computational platform. Thus, we consider here whether it is possible to capture all features of the solution, both coarse and fine, on coarse meshes. In other words we endeavor to develop a numerical scheme that shows high accuracy on crude discretizations in an effort to reduce the computational overhead associated with large scale problems that arise in higher dimensions. With this in view, let us assume that the total solution f is decomposable into f ¼ f h þ f e; where f h is the part of the solution that can be obtained with a good numerical scheme and f e is the part of the solution that is lost because it is finer than the artificial length scale h. In other words, f e represents the error in the solution. We now present a systematic procedure wherein we build this part back into the formulation. The beauty of the procedure is that it automatically yields a stabilized form of the problem. In addition, we get an explicit expression for the stabilization parameter. It is important to note that since a stabilization term is not being a priori designed for the problem, the method is amenable to approximation in n-dimensional space. Thus, it is ideally suited to Kolmogorov equations in higher dimensions, where generalization of existing stabilization methods is not straightforward and often not effective. 2. The Fokker–Planck equation The Fokker–Planck, or forward Kolmogorov, equation corresponding to an n-dimensional dynamical system subjected to Gaussian white noise excitation is given by ! n n X m X of o 1 X o ¼ ðzj f Þ þ ½H ij f ; ð1Þ ot oxj 2 i¼1 j¼1 oxi oxj j¼1 where f(x, tjx0, t0) is the transition probability density function, vector x represents an n dimensional space, and z(x) 2 Rn, H(x) 2 Q Rn · Rm are the drift vector and diffusion matrix, respectively. The initial condition is n given by f ðx; 0jx0 Þ ¼ i¼1 dðxi x0i Þ. Let def
r ¼ gradient operator in ‘n’ dimensions; def
r ¼ divergence operator in ‘n’ dimensions; def
D ¼ Laplace operator in ‘n’ dimensions:
1516
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
Furthermore, let X Rn be an open bounded region in n dimensions with piecewise smooth boundary C. We assume n P 2. The Fokker–Planck equation can be written in direct notation as f_ ¼ r ðzf Þ þ 12DðHf Þ in X 0; T ½:
ð2Þ
Eq. (2) together with the initial condition constitutes the strong form of the problem. Now, let w H1(X) be the weighting function for the probability density function. Taking the inner product of the weighting function with (2) and applying the divergence theorem yields the weak form of the problem. ðw; f_ Þ þ ðw; r ðzf ÞÞ þ 12ðrw; Hrf Þ ¼ 0 in X 0; T ½;
ð3Þ
where (Æ, Æ) is the L2 inner product. In the final step, we let V H 1 ðXÞ \ C 0 ðXÞ represent the space of piecewise continuous trial solutions and weighting functions. Substituting the discrete counterparts of f and w into (3) yields the standard Galerkin form of the problem. Remark. In general, the standard Galerkin form requires a stabilized method to yield a solution. We first present the salient features of the multi-scale decomposition that we will employ in the present development. The Variational Multiscale Method was first proposed by Hughes [13].
3. Multiscale decomposition We consider the bounded domain X discretized into numel non-overlapping regions Xe (element domains in n dimensions) with boundaries Ce, e = 1, 2, . . . , numel such that X¼
n[ umel
Xe :
e¼1
We denote the union of element interiors and element boundaries by X 0 and C 0 , respectively, X0 ¼
n[ umel
ðintÞXe
ðelement interiorsÞ;
ðintÞCe
ðelement boundariesÞ
e¼1
C0 ¼
n[ umel e¼1
and we assume an overlapping sum decomposition of the PD (probability density) field into coarse, or resolvable-scales and fine or unresolved-scales. f ðxÞ ¼ f ðxÞ þ f 0 ðxÞ : ð4Þ |ffl{zffl} |ffl{zffl} coarse scale
fine scale
Likewise, we assume an overlapping sum decomposition of the weighting function into coarse and fine and w 0 , respectively scale components, indicated by w ðxÞ þ w0 ðxÞ : wðxÞ ¼ w |ffl{zffl} |fflffl{zfflffl} coarse scale
ð5Þ
fine scale
We further assume that the fine scales, although non-zero within the elements, vanish identically over the element boundaries
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
f 0 ¼ w0 ¼ 0
on C0 :
1517
ð6Þ
We specify a direct sum decomposition on the spaces of functions for the coarse and fine scale fields V0 ; V¼V
ð7Þ
in (7) is the space of trial solutions and weighting functions for the coarse scale PD field, and is where V identified with the standard finite element space C 0 ðXÞ \ V; f 2 V
e Þ ¼ Pk ðXe Þ: VðX
Here, Pk ðXe Þ denotes the complete polynomials of order k over Xe. On the other hand, various characterizations of V0 are possible, subject to the restriction imposed by the and V0 to be linearly independent. Consequently, in the discrete stability of the formulation that requires V 0 case, V can contain various finite dimensional approximations; e.g., bubble functions or p-refinements that satisfy Eq. (6), i.e., f 0 2 V 0 ¼ ff 0 jf 0 ¼ 0 on C0 g:
ð8Þ
4. The multiscale problem We now substitute the trial solution (4) and the weighting function (5) in the standard weak form (3). This becomes the point of departure from the conventional Galerkin formulation 0 ð w þ w0 ; f_ þ f_ Þ þ ð w þ w0 ; r ðzðf þ f 0 ÞÞÞ þ 12ðrð w þ w0 Þ; Hrðf þ f 0 ÞÞ ¼ 0:
ð9Þ
Remarks (1) With suitable assumptions on the fine scale field, as stipulated in (6), and employing the linearity of the weighting function slot, we can split the problem into coarse and fine scale sub-problems, indicated as and W0 , respectively. W (2) We assume that f 0 is represented by piecewise polynomials of sufficiently high order, continuous in x but discontinuous in time. In particular, f 0 is assumed to be composed of piecewise constant-in-time functions [15]. Therefore f_ ¼ f_ . The coarse scale sub-problem, W ð w; f_ Þ þ ð w; r ðzðf þ f 0 ÞÞÞ þ 12ðr w; Hrðf þ f 0 ÞÞ ¼ 0:
ð10Þ
contains only the coarse scale components. The weighting function slot in W The fine scale sub-problem, W0 ð w; f_ Þ þ ð w; r ðzðf þ f 0 ÞÞÞ þ 12ðrw0 ; Hrðf þ f 0 ÞÞ ¼ 0:
ð11Þ
The weighting function slot in W0 contains only the fine scale components. The general idea at this point is to solve the fine scale problem to obtain the fine scale solution f 0 which can then be substituted into (10), thereby eliminating the fine scales yet retaining their effect.
1518
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
5. Solution to the multi-scale problem 5.1. Solution to the fine scale problem Consider now the fine scale part of the weak form W0 , which because of the assumption on the fine scale space, is defined over the sum of element interiors indicated as X 0 . Exploiting linearity of the solution slot on the left hand side of (11), we have 1 1 ðw0 ; f_ ÞX0 þ ðw0 ; r ðzf ÞÞX0 þ ðw0 ; r ðzf 0 ÞÞX0 þ ðrw0 ; Hrf ÞX0 þ ðrw0 ; Hrf 0 ÞX0 ¼ 0: |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 2 2
ð12Þ
Term 3
Further, consider Term 3 alone in Eq. (12) and apply the mathematical identity r ðzf 0 Þ ¼ f 0 div z þ z rf 0 ; so that ðw0 ; r ðzf 0 ÞÞX0 ¼ ðw0 ; f 0 div zÞX0 þ ðw0 ; z rf 0 ÞX0 :
ð13Þ
Now, apply integration by parts to the left-hand side of (13) ðw0 ; ðz nÞf 0 ÞC0 ðz rw0 ; f 0 Þ ¼ ðw0 ; f 0 div zÞX0 þ ðw0 ; z rf 0 ÞX0 : 0
ð14Þ 0
In (14), the first term on the left-hand side is zero because w = 0 on C (see Eq. (6)). For a piece-wise constant projection of z div z ¼ 0: Accordingly, the first term on the right-hand side is zero, and we see that the right-hand side is equal to the negative of the left-hand side. (Note: Switch positions of w 0 and f 0 to see this.) This proves that the third term in Eq. (12) is zero. Accordingly, rearranging terms in (12) yields 1 ðrw0 ; Hrf 0 ÞX0 2
¼ ðw0 ; f_ ÞX0 ðw0 ; r ðzf ÞÞX0 12 ðrw0 ; Hrf ÞX0 ¼ ðw0 ; f_ þ r ðzf Þ 12DðH f Þ ÞX0 : |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
ð15Þ
¼r
Remark. Consider the bracketed term on the right-hand side in (15) and compare it with (2). It is clear that this is the residual of the coarse scales for the Euler–Lagrange equation of the problem i.e., the fine scales are in fact driven by the residual of the coarse scales. Having now shown this, our objective is to solve (15) to get an expression for f 0 that can be substituted into (10) (the coarse scale problem), thereby eliminating f 0 yet retaining its effect. Without loss of generality, we assume that the fine scales are represented via bubble functions over the element domain; i.e., f 0 jXe ¼ be be
on Xe ;
ð16Þ
w0 jXe ¼ be ce
on Xe ;
ð17Þ
where be represents the bubble shape function over element domains, and be and ce represent the coefficients for the fine scale trial solutions and weighting functions, respectively. Substituting (16) and (17) into the fine
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1519
scale problem (15), taking the vectors of constant coefficients be and ce out of the integral expressions, and employing the arbitrariness of ce results in Z 2 e e b ¼R b r dX : ð18Þ T ðrbe Hrbe Þ dX X0 X0 We can now reconstruct the fine scale field over the element domain using (16). " Z # 2 f 0 ðxÞ ¼ be R ber dX : eT e 0 X 0 ðrb Hrb Þ dX X
ð19Þ
For multi-linear elements, the spacial part of the residual is constant. To simplify the expression we consider the constant projection of the f_ term. Accordingly, the fine scale field f 0 (x) can be written in a concise form as f 0 ðxÞ ¼ sr; where s¼R
2be X0
R
X0 eT
ð20Þ be dX
ðrb Hrbe Þ dX
ð21Þ
:
It is important to note that the matrix H possesses only one non-zero entry, i.e., the term in the hmm slot. Therefore, the final form of the stabilization parameters for the Fokker–Planck equation is s¼R
R 2be X0 be dX : ðbe;m jhmm jbe;m Þ dX X0
ð22Þ
Remarks (1) It can be proven that stabilization parameter s is a positive bounded function. (2) It is important to note that the choice of the bubble function only affects the definition of the stabilization parameter s. Given the expression for f 0 (the fine scale field), we next substitute it into (10) and get the modified form for W. 5.2. Solution to the coarse scale problem Starting with (10) and employing the linearity of the solution slot, we expand to get ð w; f_ Þ þ ðð w; r ðzf ÞÞ þ ð w; r ðzf 0 ÞÞ þ 12ðr w; Hrf Þ þ 12ðr w; Hrf 0 ÞÞ ¼ 0:
ð23Þ
Applying integration by parts to the third and fifth terms, and applying (6) gives 1 1 Þ; f 0 Þ ¼ 0: w; Hrf Þ ðDðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ ðr ðz wÞ; f 0 Þ þ ðr |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 2 2 |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} Term 3
Term 5
ð24Þ
1520
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
Combining terms 3 and 5 and substituting f 0 from (20) in (24) Þ; srÞ ¼ 0: w; Hrf Þ þ ðr ðz wÞ þ 12DðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ þ 12ðr Substituting r from (16) yields the modified formulation Þ; s f_ þ r ðzf Þ 12DðH f Þ ¼ 0: w; Hrf Þ þ r ðz wÞ þ 12DðH w ð w; f_ Þ þ ð w; r ðzf ÞÞ þ 12ðr
ð25Þ
ð26Þ
Remarks (1) The first three terms in (26) are the Galerkin terms. (2) The fourth term in (26) has appeared because of the assumption of the existence of fine scales in the problem. This term is not present in the classical Galerkin formulation for the Fokker–Planck equation, and represents f 0 in the problem. (3) The fourth term is in fact a residual-based term; i.e., the trial solution slot contains the residual of the Euler–Lagrange equations for the coarse scales. Accordingly, the formulation is consistent and accommodates the exact solution. (4) The additional (fourth) term enhances the stability of the formulation, without upsetting the consistency. It is important to note that stability and consistency lead to convergence. (5) The consistent mass matrix in this formulation gets an additional contribution from the stabilization term Þ; sf_ : ð w; f_ Þ þ r ðz wÞ þ 12DðH w ð27Þ (6) For the case of linear elements, a comparison of the modified formulation with the two celebrated staÞ bilized methods, i.e., the SUPG and the GLS methods is straightforward. If we drop the þ 12 DðH w term from the weighting function slot of the stabilization term, we recover the SUPG formulation Þ term in the weighting function for the Fokker–Planck equation. If we change the sign of þ 12 DðH w slot of the additional stabilization term and add the weighting counterpart for the rate term, we recover the GLS formulation for the Fokker–Planck equation. However, it is important to note that the stabilization terms in SUPG and GLS contain stabilization parameters that designed based on the consideration of achieving nodally exact solution in the one-dimensional setting. Generalization of the nodally exact stabilization parameters in 1-D to the case of nodally exact parameters in multi-dimensions has been a challenging task. A literature review in the field of fluid mechanics would reveal various attempts in designing the optimal stabilization parameters. (7) In the proposed formulation the definition of the stabilization parameter s appears naturally and is given in Eq. (21). This is in sharp contrast to the SUPG or the GLS stabilization schemes for the Kolmogorov equations. (8) It is important to note that the use of bubble functions resides completely in the definition of s, as given in Eq. (22). (9) The modified formulation is completely expressed in terms of the coarse scale PD function f . This means that the dimension of the element stiffness matrix is the same as that for the Galerkin method. (10) The proposed formulation is stable in Rn, where n is the spatial dimension of the Fokker–Planck equation considered. (11) The formulation accommodates Lagrange-based C0 functions. In compact notation, the multi-linear shape function can be expressed as
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
N a ðnÞ ¼
n 1 Y ð1 þ nai ni Þ 2n i¼1
1521
ð28Þ
n ¼ n1 ; n2 ; . . . nn : (12) The formulation also accommodates triangles, tetrahedron and their generalizations in n-dimensions.
6. Numerical examples We have developed 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals (see Fig. 1) for application to the Fokker–Planck equation. We have employed the simplest polynomial bubbles, which in the element natural coordinate frame are as follows: be ðn; gÞquadrilateral ¼ ð1 n2 Þð1 g2 Þ;
ð29Þ
be ðn; gÞtriangle ¼ 27ngð1 n gÞ:
ð30Þ
6.1. The rate of convergence study The first numerical simulation is a study of convergence rates on a model problem. The domain under consideration is a bi-unit square with origin at (x, y) = (0, 0). The exact solution is assumed to be of the form f ¼ sin
2px 2py sin : L L
ð31Þ
The ‘‘forcing function’’ as evaluated from Eq. (1) is given by 2p 2px 2py 2px 2py p 2px 2py z1 cos sin z2 sin cos sin sin f_ ¼ : L L L L L L L L
Fig. 1. A family of 2-D linear and quadratic elements.
ð32Þ
1522
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
In specifying the boundary-value problem, the forcing function is integrated over X while f = 0 is prescribed at the boundary. Note that H22 = 1 in this convergence rate study. The meshes employed in the convergence study are regular structured meshes. For linear quadrilateral elements, the meshes consist of 100, 400, 1600 and 6400 elements. The linear triangular element meshes consist of exactly twice as many elements and are constructed by bisecting the quadrilaterals along the main diagonal. For quadratic quadrilateral elements, the meshes employed consist of 100, 400, 1600 and 6400 elements. Again, the quadratic triangular element meshes consist of twice as many elements. The element mesh parameter, h, is taken to be the edge length of the elements for the quadrilaterals, and the short-edge length for triangles. Convergence rates for jzij = 10 are presented in Figs. 2 and 3. Fig. 2 shows the L2-norm and H1-seminorm of the scalar field for the linear quadrilateral and triangle. We see that optimal rates are achieved. Fig. 3 presents the corresponding rates for the quadratic elements. Once again we see optimal convergence rates. Fig. 4a and b present the elevation plots of the exact and the computed fields for a representative
Fig. 2. Convergence rates for linear elements.
Fig. 3. Convergence rates for quadratic elements.
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1523
Fig. 4. (a) Elevation plot of the scalar field (exact solution). (b) Elevation plot of the scalar field for the 4-node quadrilaterals (20 · 20 element mesh).
mesh of 4-node elements. These plots present a qualitative comparison of the computed solution with the corresponding exact solution, and show a stable response of the elements. 6.2. The linear oscillator (2D calculations) Simple verification of the proposed formulation (26) is afforded by the solution of the Fokker–Planck equation for the two-dimensional linear oscillator subjected to additive white noise excitation, for which an exact solution exists. The non-dimensionalized Fokker–Planck equation is given by 1 of o2 f o o ¼ 2n 2 ð2ny þ xÞf ðyf Þ: 2p os oy oy ox
ð33Þ
Here, x = x1/r, y = x2/rx0, s = x0t/2p, and r2 ¼ pS 0 =2nx30 . The damping parameter is n = 0.2. The computational domain is given by 10 6 x, y 6 10; the extent of the domain, 10, is chosen based on experience as
0.16 0.14
fx (0,0,τ)
0.12 0.10 0.08 0.06 0.04 0.02 Present Calculations 0.00 0
2
4
6
8
10
12
14
16
18
time Fig. 5. Probability density function for the linear oscillator at (0, 0).
1524
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
necessary in transient computations with non-zero initial conditions to minimize the incidence of probability mass reflecting off the boundaries and back into the interior. The mesh employed in the simulation is composed of 50 · 50 4-node elements, and the time-integration method is the Crank–Nicholson algorithm with Dt = 0.005. Boundary conditions are set such that f = 0 at the boundaries. The initial condition is a Gaussian distribution with standard deviation of unity, centered at (x, y) = (3, 3). Fig. 5 shows the evolution of the probability density function at (0, 0). Note that, at stationarity, the value of the distribution at the origin converges to f(0, 0) = 1/2p, the exact solution. Fig. 6 shows the elevation plot of the probability density function at times 1, 2 and 4 s, superimposed on common axes. Fig. 7 shows the exact probability density function at stationarity, while Fig. 8 shows the corresponding computational result. There is no observable difference between the two to the machine precision.
Fig. 6. Evolution of the probability density function for the linear oscillator.
Fig. 7. Exact probability density function for the linear oscillator.
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
1525
Fig. 8. Computed converged probability density function for the linear oscillator.
7. Conclusions We have developed a multi-scale formulation for the Fokker–Planck equation of a class of oscillators in Rn. The key point of the method is a multi-scale decomposition of the field variable into resolvable/coarse and unresolvable/fine scales. The proposed multi-scale method leads to a stabilized formulation which is free of any user-defined parameters. Since the fine or the subgrid scales in the formulation are proportional to the grid scale residual, the method is consistent. We have tested the formulation on linear and quadratic triangular and quadrilateral elements and have obtained optimum numerical convergence rates in the norms considered. The formulation was applied to a model problem, the two-dimensional linear oscillator subjected to external white noise excitation, for which an exact solution exists. The computational result is a converged solution on a mesh of 2500 four node quadrilateral elements, representing an order of magnitude savings over previously reported converged results. While the two-dimensional result is not overly compelling, it does show significant promise for higher dimensional calculations. Prior FE solutions required on the order of 100 elements per spatial dimension for nominal convergence of a well-behaved dynamical system. Thus, a four-dimension problem which typically required 100 million elements would, incorporating the multiscale formulation, require nearly two orders of magnitude fewer. For the moment, efforts continue toward further refinement of the procedure with eventual application to larger scale systems subjected to both additive and multiplicative excitations.
References [1] [2] [3] [4] [5]
G.I. Schueller (Ed.), A state-of-the-art report on computational stochastic mechanics, Prob. Engrg. Mech. 13/4 (1997) 197–321. J.D. Atkinson, Eigenfunction expansions for randomly excited non-linear systems, J. Sound Vib. 30 (2) (1973) 153–172. Y.K. Wen, Approximate method for nonlinear random vibration, J. Engrg. Mech. Div., ASCE 101 (EM4) (1975) 389–401. Y.K. Wen, Method for random vibration of hysteretic systems, J. Engrg. Mech. Div., ASCE 102 (EM2) (1976) 249–263. L.A. Bergman, J.C. Heinrich, On the reliability of the linear oscillator and systems of coupled oscillators, Int. J. Numer. Methods Engrg. 18 (1982) 1271–1295. [6] L.A. Bergman, B.F. Spencer Jr., Robust numerical solution of the transient Fokker–Planck equation for nonlinear dynamical systems, in: N. Bellomo, F. Casciati (Eds.), Nonlinear Stochastic Mechanics, Springer-Verlag, Berlin, 1992, pp. 49–60. [7] B.F. Spencer, L.A. Bergman, Numerical solutions of the Fokker–Planck equations for nonlinear stochastic systems, Nonlinear Dyn. 4 (1993) 357–372.
1526
A. Masud, L.A. Bergman / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526
[8] S.F. Wojtkiewicz, L.A. Bergman, B.F. Spencer, Numerical solution of some three-state random vibration problems, in: S.C. Sinha (Ed.), Vibration of Nonlinear, Random, and Time-Varying Systems, ASME DE-84-1, 1995, pp. 939–947. [9] W.V. Wedig, Utz Von Wagner, On the calculation of stationary solutions of multi-dimensional Fokker–Planck equations by orthogonal functions, Nonlinear Dyn. 21 (2000) 289–306. [10] T.J.R. Hughes, A.N. Brooks, Streamline-upwind/Petrov–Galerkin methods for advection dominated flows, in: Proceedings of the Third International Conference on Finite Element Methods in Fluid Flow, Banff, June 1980, pp. 283–292. [11] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1–3) (1982) 199–259. [12] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (2) (1989) 173–189. [13] T.J.R. Hughes, Multiscale phenomena: Greens functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. [14] A. Masud, T.J.R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg. 191 (39–40) (2002) 4341–4370. [15] A. Masud, On a stabilized finite element formulation for incompressible Navier–Stokes equations, in: Proceedings of the Fourth US–Japan Conference on Computational Fluid Dynamics, Tokyo, Japan, May 2002. [16] M. Ayub, A. Masud, A new stabilized formulation for convective–diffusive heat transfer, Numer. Heat Transfer, Part B 44 (2003) 1–23. [17] A. Masud, R. Khurram, A multiscale/stabilized finite element method for the advection–diffusion equation, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018.
