- Email: [email protected]
Petrov-Galerkin finite element solution for the first passage probability and moments of first passage time of the randomly accelerated free particle
COMPLrI'ER M E T H O D S IN A P P L I E D M E C H A N I C S A N D E N G I N E E R I N G 27 (1981) 345-362 NORTH-HOLLAND PUBLISHING COMPANY
PETROV-.GALERFAN FINITE ELEMENT SOLUTION FOR THE FIRST PASSAGE PROBABILITY AND MOMENTS OF FIRST PASSAGE TIME OF THE RANDOMLY ACCELERATED FREE PARTICLE L.A. BERGMAN Department of Theoretical and Applied Mechanics. University of Illinois. Urbana. IL 61801. U.S.A.
J.C. HEINRICH Department of Aerospace and Mechanical Engineering. The University of Arizona. Tucson. A Z 85721. U.S.A.
Received I() August 19~() R e v i ~ d manuscript received 1~ January 1981
A numerical solution to an initial boundary value problem governing the probability of failure of a randomly accelerated free particle is obtained using a Petrov-Gaderkin tinite element method. This direct solution is the first successful one. ,'rod no o~hers ha~.e bccn reported ia the literature. A solution of the Pontriagin-Vit.t equation for the t|mc to first passage of the particle is obtained first: in this case an analytical solution is available and used to e~aluate the numerical algorithm. Extensions to the solution of other stochastic differential equations, in p~:rt|cular th('Jse governing the probability of failure of the linear oscillator, and applications to ,,tructur;d dsnam~cs are discussed
I. Introduction
In recent years a number of probability based codes have been formulated as a ba,,~is for safety analysis and design, providing a consistent framework for the swtematic as,~essment of the reliability of a structure. This evolution has occurred as a result of high risk technological development in structural systems. The determination of the probability of failure of such a system during its expected life is. in the most general sense, the reliability problem. When failure is defined as the first excursion of a system beyond some limit, whether displacement, strain or stress, due to a random excitation, the problem is characteristically referred to as a first passage problem. The simplest problem in structural reliability is represented by the linear single degree-offreedom oscillator subjected to stationary Gaussian white noise excitation, for which the displacement and velocity responses form the components of a two dimensional Markov vector process [1]. Therefore, the evolution of the transition probability density function is governed by a Fokker-Planck equation and, alternately, by its adjoint the backward Kolmogorov equation. Although the latter leads to a well-posed initial-boundary value problem for the system probability of failure allowing the use of direct numerical integration procedures, prior efforts 0045-7825/81/0000--0000/$02.50 © 1981 North-Holland
346
L.A. Bergman, 2.C. Heinrich, Galerkin and randomly accelerated particles
in this direction have been unsuccessful [2]. This has led to the formulation of a variety of approximations which are discussed in [1, 3, 4] and references therein. In the present work, we study a limiting case of the linear oscillator wherein the spring and dashpot are removed leaving a randomly accelerated free particle. The resulting initial boundary value problem contains all of the numerical difficulties of the more general problem with one fewer parameter. Furthermore, an analytical solution exists [5] for the PontriaginVitt equation governing the first moment of time to first passage of the free particle, permitting an assessment of the accuracy of the algorithm. A numerical solution of the transient equation for the free particle is obtained using a Petrov-Galerkin finite element method. This solution is, to the authors' knowledge, the first of its kind. The results presented here show that the method constitutes a powerful tool for the solution of a variety of problems arising in stochastic processes for which solutions have not before been available.
2. Problem description We consider a randomly accelerated particle which satisfies the stochastic differential equation
J~(t) = n(t) t > 0 ,
(1)
where X(t) is a random displacement process, n is a stationary Gaussian white excitation, the dot represents differentiation with respect to time, and with initial position and velocity X(0) = x, X(0) = y, over a 'safe' spatial domain O = {(x, y)/lxl- 1, lyl < o0}, which is shown in fig. 1. Let P(tlx, y) be the probability of failure of the particle to remain in the safe region, having started at (x, y), and without having previously failed and re-entered the safe domain. The function P then satisfies the initial boundary value problem c3P c3P c~2p Ot = y - ~ + ¢rSo 0y2 ,
(2)
with initial and boundary conditions
POI x, y) = 0
for (x, y) E n ,
P(t I-1, y)= 1
for y < 0 ,
P(t[ 1, y ) - 1 P(tl x, y ) ~ 1
for y > O,
(3)
(4)
as lyl--' oo,
where So is the magnitude of the spectral density of the excitation. A characteristic of this typ~ of equation is that no boundary condition is prescribed on part of the boundary, in. this c~,~ along the boundaries x = - 1 , y > 0 and x = 1, y < 0. A short account of the derivation of eqs
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
347
T=O p=l
. _.Y
T=O
p=t
(-I,O)
.,
(I,0)
(o, o)
x
sT p } FREE
T--O p=l
I I
i
-y T=O. p.'I
Fig. 1. Spatial domain and boundary conditions for cqs. (2) anti (¢i).
(2)-(4) is given in Appendix A. For further details the reader is referred to [1, 2.6]. it has been shown [2] that eqs. (2)-(4)constitute a well-posed problem: however to the author's knowledge no direct solution either analytical or numerical has been reported. Defining the function T(x, y) to be the expected value of the time to first passage of the randomly accelerated particle (i.e., the expected value for the first time the particle reaches a position IX(t)l = 1), T(x, y) is then the solution of the boundary value problem (see Appendix
A) 02T
OT
with boundary conditions
T(-l,y)=O
fory
T(I, y)= 0
for y > O,
Ttx, y)--, 0
as ly[--, 00.
(¢,)
Eq. (5) is known as the Pontriagin-Vitt equation for the first moment of time to first passage of the randomly accelerated particle. Problem (5)-(6) was solved analytically and numerically by Franklin and Roclemich [51 for
348
L.A. Bergman, J.C. Heinrich, Oalerkin and randomly acceleratedparticles
the case ¢rSo= 0.5. In their numerical calculations they used two finite difference methods previously developed in [7, 8]. However, their numerical results were not satisfactory, and they only obtained reasonable accuracy by averaging the nodal solutions in one of their methods. The analyst faces two basic difficulties when attempting the numerical solution of prob)em (5)-(6). The first stems from the character of the operator, Eqs. (2) and (5) fall into the general category of degenerate-elliptic or elliptic-parabolic operators which are positive semidefinite. In fact, eq. (5) has the form of a heat equation with time x and temperature T propagating from the right when y > 0 and from the left when y < 0. This requires that the numerical scheme be capable of taking the direction of the initial velocity y into account and producing discretizations that are biased accordingly. The oecond difficulty comes from the form of the solution of problem (5)-(6). Franklin and Rodemich [5] obtained an analytical solution that along y = 0 takes the form
T(x, 0 ) =
x2)~'6[F(-1/3,1, 7/6,½(1- x))+ F(-1/3, 1, 7/6, ½(1+ x))], C--3~"6/(2F(1/3)), F is the hypergeometric function and F the gamma C(1-
(7)
where function [9]. Clearly, the first partial derivatives of T with respect to x are unbounded at x = +_1, which causes additional computational difficulties and raises questions about convergence of the numerical algorithms. These will be discussed in a later section. In what follows we will assume ~rS0= 0.5, define a Petrov-Galerkin finite element method for the numerical solution of problem (5), (6), analyze its behavior and accuracy by comparison with the analytical solution (7) along the x-axis and discuss convergence. We then extend the algorithm to the transient equation (2), which is solved, to our knowledge, for the first time and evaluate the results.
3. Numerical solution of the Pontriagin-Vitt equation The partial differential equations (2) and (5) can be classified in general as convective diffusion equations in a highly anisotropic diffusive media. The difficulties related to the numerical solution of such equations when the first derivative (convective) terms are dominant led to the development of Petrov-Galerkin finite element methods by Heinrich et al. [10, 11] which are used in the context of the present problem. For reasons to be discussed in the next section, we restrict ourselves to a bounded rectangle ~ r = {(x, y)I Ixl-< 1, lYl-< Y} where Y is a large number (see fig. 1). The asymptotic boundary condition is prescribed at lYl = Y- Following a standard weighted residuals formulation [12], we replace the problem of finding the solution T(x, y) of the differential equation (5) over ~2~, by that of finding a function T(x, y) such that
fa Y{- ~IOWOT ~ ~yy 4 - y W ~OT} - d x d y - - fa
Wdxdy.
(8)
Y
T(x, y) is then
approximated by T(x, y)----Zk Nk(x, y)Tk, where Nk(x, y) are the well known bilinear quadrilateral basis functions and Tk the approximation to T at node k. For a particular element in isoparametric coordinates (~, ,/) the basis functions are given by
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
N'(~:, r/)= ¼(1-
349
~e)(1 - I,/)
N2(se, n)= I(1 + ~e)(1 - rl) N3(se, ~)= ¼(1 + ~e)(1+ n)
-1<_~, r/_< 1.
(9)
N'(~, n) = I(1 - O 0 + n) Following Heinrich et al. [10], the weighting functions are chosen to be
W2(~, 7/)= N2(g, 7/)(1 + 3 a ( 1 - ~:)),
W'(se, n) = N'(~:, "r/)(1 --~a(1 + ~)), W3(6 r/)= N3(se, r/)(1 +-~a(] - se)),
W4(~, 1~) -- N 4(~:, 7/)(1- ~a(1 + ~:)) ,
(10)
where t~ = 1 if y < 0, 0 if y = 0 and - 1 if y > 0. Fig. 2 shows the element notation together with two of the above weighting functions for the case t~ = - 1 , that is, for a node in the upper half plane. The other two are mirror images of the ones shown about the line ~, = 1. For t~ = 1 the weights are mirror images about the line ~ = - 1 . The effect of the special weighting functions in this numerical algorithm is better understood if we examine the difference equations generated for each node. We consider a rectangular mesh with nodes (xi, yj) at which T(x, yj)= T~. For simplicity we also approximate the coefficient y of the x-derivative by yj in the equation for node (x, y,). Under these conditions we obtain: (a) For y > 0 and a = - 1:
+
+ ~
+(-2Ax-2ydty2)T~+ - T + 2 y ~ y 2 +
y,A~22)T~-.~+\4 !
T~.,+ -
+Y/ 2 ] ~'~=-3Ax'~y2" N4
w4 4
Nil WI
T',:'~+Ax-y,
T',*' (II)
N3 w3
2NZ - g W2 '
I/*
.~
\ a=-I
350
L.A.
(b) For y
=
Bersman, ZC. Hein6ch, Galerkin and randomly accelerated particles
0 and a
=
O:
(12)
+ (~)T~+-I + (Ax)T~+' + ( ~ ) T { $ 1 = - 3 A x A y 2 •
(c) For y < 0 and a = 1 (here yj = [yj]):
* (--2Ax-- 2yjAy2)T~"['I-~)T~+I'Jt'(5~X3c--~J~2~erJ+'21i-1"~" ( Ax-y]~--~-2 )T,'+I + (-- T AX~'r'/+' ! " ÷ ' = -3AxAy 2
(13)
In fig. 3a we further illustrate the nine node stencil for the case y > 0 and a = - 1. It can be seen that the Petrov-Galerkin finite element method not only approximates the first derivative term in a backward way with respect to the velocity direction bu~ it also shifts the weighting of the second order term in the x-direction in the nine node stencil, still maintaining the symmetry of the approximation in the y-direction. This property of Petrov--Galerkin schemes, to date not well understood, appears to be a key factor in the algorithm's success. Furthermore, in the presence of
Ax-yj AY22
O
O
(Xt_l,Yj+ 1)
(xi, Y:]+I)
(x.I+I,YJ+ 1)
~x --
-2~x-2y ~y2 Oj
_ 5~x4+ 2yjAy2
- 7~x4-2YJAy2
O (xi+I, Yj)
(x i ,
(Xt_l,Yj) ~x ---~ O (xi_ 1, Yj_l)
(xt, Yj) 4y2 ~x-yj 2 O (xi, Yj_I)
a
5Ax8+ YJAy22
7Ax _Ay2 s "Yl 2
- -'8Az
O
8
YJ 2 O
(xi+1,
Yj_I)
5Ax8+ YJAy22
O'
0
(xl] YJ+I)
(xi+1,
/j+l)
_ 5~x4j 2YJAyz
O
O (xi+ 1, Yj)
Yj)
7Ax Ay2 8 - YJ 2 O (xi, Yj_l)
5Ax [ 8
~y2
YJ 2
O
(xi+I, Yj_I) b
Fig. 3. Differences stencils for Petrov-Galerkin Finite Element Method; (a) full molecule for interior node, (b) molecule for a node at the free left-hand side boundary.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
351
non-constant source terms, as is the case if higher moments are to be computed using the Pontriagin-Vitt equation [13, 4], it can be easily proved that the limiting one-dimensional case dT/dx = Q,
(14)
T(O)= 0
is exactly integrated when O is linear [14]. Any other existing backward differencing scheme applied to equation (14) will give a numerically diffused solution in the presence of a linear source term O, the difference in the Petrov-Galerkin scheme being that O is weighted consistently, with the consequent re-distribution of the nodal contributions of the source term to the right hand side vector. In fig. 3b we show the contributions to the stiffness matrix for the case of a node in the part of the boundary at x = - 1 where no boundary condition is imposed. It is clear, since the derivative of T with respect to x is approximated by a fully backward difference with respect to the velocity, that the numerical scheme does not impose any kind of boundary condition there. In fact, the present scheme is free of any kind of natural boundary condition. The resulting discretization is assembled in the form (15)
KT= B,
where
k° =
v
1 8W' aN__~ j + yW~
2 8y
c~y
Ox J
dx dy
and
b,
~,
dx d y .
This system of linear equations is solved by a direct frontal elimination procedure for unsymmetric matrices [15].
4. Convergence of the Petrov--Galerkin method Several questions arise concerning convergence of the numerical algorithm presented in the previous section. At this point, it is not within our capabilities to present a proof for the convergence of the method. In fact, none of the, existing theorems appears to be applicable in this case. Our arguments will be later substantiated by the numerical evidence. The fact that no second derivative term with respect to x appears in eq. (5) implies that c~ does not go to zero as the mesh size decreases. This makes the results of Griffiths and Lorenz [16] not applicable to our equations. A straightforward truncation error analysis applied to eqs. (11)-(13) indicates that the truncation error is O ( A x + dy:) as expected. However the leading terms depend on the higher order derivatives ~3Y/o~x~y 2 ,
02T/~x 2
and
8"T/Sy 4 .
Examination of equation (7) reveals that along the x-axis the solution T is the product ~,f an even function of x by ( l - x 2 ) TM, and the first partial derivatives with respect to x are unbounded in a neighborhood of x = _+1, y = 0. Furthermore the square of the first derivative
352
L.A. Bergman, .I.C. Heinrich, Galerkin and randomly accelerated particles
is non-integrable. As a consequence, in this particular case the truncation error in (11), (13) is unbounded and from the Petrov-Galerkin point of view we are violating the basic requirement of the solution lying in the space of functions with a finite energy. This was resolved by Franklin and Rodemich [5] by constructing a family of similar functions T(x, y, a) such that (a) T(x, y, a) is smooth for Ix[--< 1, [y[-< Y, if a > 1, (b) T(x, y, a)--, T(x, y) as a -- 1, and (c) LT(x, y, a ) ~ LT(x, y) as a -~ 1, where L =-y(O/Ox)-~(0Z/0y:). These functions essentially shift the singularity to the points x = +-a, y = 0 with ial > 1 as a consequence though no convergence rates will be obtainable. For more details see Franklin and Rodemich [5]. The fact that the domain is unbounded in the y-direction presents no difficulty from the finite element point of view since asymptotically the solution has the behavior 1-x Y
T(x, y)
as y-,0o,
T(x, y ) ~ - l~+ x Y
as y - , - 0 o ,
(16)
and this is square integrable over the infinite region. However, in order to avoid dealing with an infinite region and unnecessary complications introduced by the use of 'infinite elements' (see e.g. [12]) we choose to truncate the region at a far enough distance Y and accept introducing yet another small error in the manner of Franklin and Rodemich [5]. It will be clear later from the numerical results that the error introduced by this approximation is extremely localized and has no significant effect even if homogeneous boundary conditions are imposed at the truncated boundary for fairly small values of Y.
5. Numerical solution of the transient equation
We now turn to the solution of problem (2)-(4). The spatial operator remains unmodified from the previous case so now the Petrov-Galerkin finite element method is applied to eq. (5) with the right hand side replaced by the time derivative. This yields the semi-discretized system M//-KP
= 0,
(17)
where the stiffness matrix K is the same as in eq. (15), P is the vector of nodal approximations Pk to P(tlxE, Yk) and/~ is its first time derivative. The mass matrix M is now given by
m° = fn WiNS dx dy.
(18)
Y
Following the notation of section 3, the resulting assembled mass matrix for a typical node (xi, Ys) is shown in fig. 4a. We point out that to make it compatible with the expressions for the stiffness matrix of fig. 3 the terms in the mass matrix must be multiplied by 3Ay. An interesting situation arises when a boundary node with no boundary condition applied is
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles AxAyr~ a 12 (3 + ~)
AxAy.I
AxAy
O'
Yj+l)
9 O
~C i 12 "3 + 0
a 2)
4AxAy
I
AxAy.2
-12 ( 3 '0
0
(xi_I, Yj_I)
(x t , Yj_I)
I,
AxAy:! a O
0
(xi+I, Yj)
AxAy.1
(Xi+l~ Yj+I)
3 t~ - ~)
o
AxAv 9
O
[ Y]+I)
(xi'
(xi[ YJ)
A:'~y.l a
ct
O
O
(xi+I, Yj+I)
9 O
(xi_1, Yj)
AxAy.2
(xi, Yj+!)
(xi-Ii
nxAy.t a 3 ~3 + ~) O
a
n 2)
(xi+1, Yj_!)
a
(xi+I, Yj)
(xii YJ)
AxAy.2 .,~ 12 (3 - 2)
AxAy I - ~) 12 (3
0 (x.i ,
353
0
Yj_!)
(xi+I, Yj_I) b
Fig. 4. Mass matrix for Petrov-Galerkin Finite Element Method; (a) interior node. (b) node at the free left-hand side boundary.
reached, as shown in fig. 4b for the case y > 0. Since a = - 1 in this case, the method will increase the mass at these boundaries. This should result in a stabilizing effect there, overall conservation being achieved by subtraction of the same amount at the opposite boundary. This, however, indicates that a mesh symmetric about the y-axis should be used. It is clear in fig. 4 that the Petrov-Galerkin method redistributes the masses of the non-central nodes to the point of assigning negative masses to those downstream with respect to the velocity direction, while keeping the symmetry in the y-direction. Eq. (17)is further discretized using the theta method [12]: that is [M-
At(1 - 0)Kle
"+' = [ M -
AtOK]P'.
(l.)
where 0 <-0-< 1, n denotes the current time and n + 1 the time to advance to. For 0 = 0.5 the second order in time Crank-Nicolson marching scheme is obtained, and for 0 = 0 the fully implicit forward marching method, both unconditionally stable. These two methods have been used for the time discretization with time steps varying in the range 0 . 0 ~, ~ A t ~ O.2 s. T h e critical time step, beyond which a fully explicit (0 = 1.0) .scheme would be unstabie, was estimated in accordance with Myers [20] to be 0.025 s. An upper bound for the time steps prior to excessive loss of accuracy or the onset of oscillations was found to be dependent on the problem and in particular on the mesh.
6. Numerical results for the Pontriagin-Vitt equation Numerical solutions to problem (5), (6) have been obtained fol a var:~ety of situations to assess the effect of the mesh, boundary conditions and the artificial boundary y = ± Y.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
354
The results indicate that the boundary conditions at y = - Y have a very localized effect, since we are interested mainly in the results about y = 0 where the maximum gradients occur. We find that for a value of Y as small as 5 no difference can be detected in the solution along y = 0 when the asymptotic boundary conditions (16) are replaced by homogeneous conditions. A second question is how the imposition of a homogeneous boundary condition at the points (+-1, 0) affects the solution (Franklin and Rodemich [5] did not impose T(+I, 0)= 0 in one of their finite differences schemes). It was found that although the imposition of T(+__l, 0)= 0 causes some oscillation close to those nodes, this oscillation disappears rapidly (within two nodes) and can be reduced considerably by adequate mesh refinement in the neighborhood of the singular points. The relaxation of this condition on the other hand results in a small overall loss of accuracy, which is minimized through mesh refinement, however no oscillations are observed in this case. From the above it also becomes clear that the best choice of mesh is an irregular one refined in the neighborhood of the points (x, y ) = ( ± l , 0), regardless of the boundary condition strategy. Here we present the results obtained with a 20 by 20 element mesh, with nodes symmetrically distributed about the x and y axis with coordinates at x = 0.0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.975, 1.0 and at y = 0.0, 0.2, 0.4, 0.7, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0. The artificial boundary is at Y = 5, and homogeneous boundary conditions were imposed there. In table 1 we show results along the x-axis obtained with and without imposing the conditions T(+_l, 0)= 0. These are compared with the analytic solution and with the finite difference solutions of Franklin and Rodemich [5] using an 11 by 21 nodes mesh with Y = 20. Fig. 5 shows the finite element solution obtained imposing the condition T(_+I, 0)= 0, for initial velocities y = 0.0, 0.4, 1.0 and 2.0; and in fig. 6 we show the three dimensional surface for the entire solution. The agreement of the finite element calculations with the analytical solution is excellent (exact to three significant figures) as shown in table 1. We point out again that refinement about the points (x, y) = (±1, 0) plays a critical role in the accuracy of the solution, and nothing can be said about convergence rates there. Table 1 Analytical and numerical results for the Pontriagin-Vitt equation X
EXACT
FEI a
FE2 b
FD1 c
FD2 ~
X
EXACT
FD3 "
0.0 0.2 0.4 0.6 0.8 0.9 1.0
2.312 2.289 2.216 2.076 1.813 1.580 0.000
2.314 2.290 2.217 2.077 1.826 1.578 0.000
2.317 2.29g 2.220 2.080 1.828 1.598 0.762
2.151 2.125 2.042 1.878 1.554 -0.000
2.475 2.115 2.409 1.850 2.171 m
0.000 0.182 0,364 0.545 0.727 0.909
2.312 2.293 2.233 2.123 1.933 1.550
2.319 2.268 2.246 2.088 1.980 1.432
"FE1 = Finite elements with T ( - 1 , 0) = 0. Finite elements with T(4-1, 0) free. CFD! = Biased finite differences. dFD2 = Biflanced finite differences (no averaging). t "FD3 = Balanced finite differences (averaged). bFE2 =
Ref. [5].
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
355
~ " , , ~ y =0.0 2.0 I..-
O0 INITIAL DISPLACEMENT x
Fig. 5. Expected values of time to first passage T vs. initial displacements for initial velocities y = O.O, 0.4. 1.0 and 2.0.
x
o
•
v
,
•
~
Fig. 6. Three dimensional surface for the expected value of time to first passage T
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
356
7. N u m e r i c a l r e s u l t s f o r the transient e q u a t i o n
Considering the previous experience in the solution of the Pontriagin-Vitt equation (5), numerical experiments for the solution of the initial boundary value problem (2) were restricted to the case Y = 5. We present results of computations performed in a 196 element irregular mesh with nodes symmetrically distributed about the x and y axis and coordinates Ix[ = 0.0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1.0 and [y[ = 0.0, 1.0, 2.0, 3.0, 3.5, 4.0, 4.5, 5.0. Table 2 shows results obtained using both the Crank-Nicolson and fully implicit methods, for the cases P(t[_+1,0)= 1 and P(ti -_.1, 0) free. In the same table we show some values for T(0, 0) obtained by calculation of the first moment of the numerical solution i.e.
Ktg(tlx, y;n)dt,
T(O,O)=
(20)
w h e r e g ( t l x , y ; D ) is t h e p r o b a b i l i t y density of t i m e to first passage (see A p p e n d i x A). C o m p a r i s o n with t h e exact ~Jiation s h o w s less t h a n 5 % error. Fig, 7 s h o w s t h e probability of failure at t h e p o i n t (0, 0) versus time, a n d it is c o m p a r e d with t h e results o b t a i n e d using a t w o m o m e n t m a x i m u m e n t r o p y distribution [17] with o r d i n a r y m o m e n t s nl = 2.295 a n d n2 = 7,29, a n d an e x p o n e n t i a l (one p a r a m e t e r ) distribution [6, 13]. A n excellent overall a g r e e m e n t can b e
Table 2 Evolution of the probability of failure at the origin P(z [0, 0); Petrow-Galerkin finite element solution
P(t tO, O) t (s) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
1"
2b
3c
0 . 0 0 0 0.000 0 . 0 8 9 0.062 0.462 0.484 0.707 0.722 0 . 8 3 7 0.848 0 . 9 1 0 0.917 0 , 9 5 0 0.955 0 . 9 7 3 0.97.6 0 . 9 8 5 0.987 0.992 0.993 0 . 9 9 5 0.996
T(0, 0) [Exact = 2.312]
"At = 0.2, 6' -- 0.0 bat = 0.1, 0 = 0.0 ~At = 0.05. 0 = 0.0
0.000 0.042 0.497 0.729 0.853 0.921 0.958 0.977 0.988 0.993 0.996
A.d
0.000 0.000 0.495 0.728 0.85J, 0,9"23 0.959 0.978 0.989 0.994 0.997.
2.483
P(t [ +--1,0) Free.
dAt = 0-1, 0 = 0-5 J
~At = 0.1, O = 0.0 / fat = O.1, 0 = 0.5 f
P(t l +-1,O)= 1.
5e
0.000 0.065 0.506 0.748 0.870 0,933 0.966 0.983 0.991 0.995 0.998
6f
0.000 0.000 0.517 0.756 0.877 0.939 0.969 0.985 0.992 0.996 0.998 2.405
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles m
1.0
..-- -"
f
//
~
""
" -
. . . . . . . . . . . . . . . .
357 , ° .
o°'°"
.....-""
°°°°°
0.8
...~r/
n~ _J
-
0.6
..""/y •"
h
0 )-
_J ,~ .ff
ID 0Iz
..."
0.4
-
/
77
.""
FEM
0.2
O0 0
2
4
6
8
~0
TIME (SEC)
Fig. 7. Probability of failure P vs. time for initial conditions of zero displacement and velocity.
IO
- -
08 ;
iAL DISPLACEMENT AND vELOCF~;Y
06
(~,y} = ( - 0 6 , 0
.. " 04
-
:(-04,
00,~
/
"
I
;
/
.
I
."
(=,y}
/
.
0,2
.
O)
/ //
oo
o
2
4
6
8
tO
T I M E (SEC)
Fig. 8. Probability of failure P vs. time for three different nonhomogeneous initial conditions of displacement and velocity.
358
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
observed in fig. 7 with a significant correction to the approximations shown by the finite elemev.t solution for lower times 0 < t < 1.5. In fig. 8 similar curves are shown for the points (0, 2), (-0.4, 0), (-0.6, 0). These plots come from run Number 3 in table 2. Fig. 9 shows the evolution of the three dimensional surface with time. After seven seconds the change in the value P(t 10, 0) become less than 0.2% every 0.1 second, and changes are no longer observable. In general few difficulties were encountered in the solution of the present problem. The initial time steps were subject to some oscillatory behavior due to the singularity in the boundary conditions. The fully implicit method (0 = 0) shows less disturbances than the Crank-Nicolson method (0 = 0.5) as would be expected, although the latter yields what appears to be a more accurate solution at later times. Also, in table 2 we observe differences
o o,>c, ~''~'~
a
t
=0 l se¢
t = 0 0 sec
c
t :05sec
t --I O s e c
t= t=2Osec.
f = 4.0see
t ='/'.0 $eC.
Fig. 9. T i m e e v o l u t i o n of t h e t h r e e d i m e n s i o n a l p r o b a b i l i t y of f a i l u r e surface; (a) initial c o n d i t i o n s , (b) t = 0.1 s, (c) t = 0.5 s, (d) t = 1.0 s, (e) t = 1.5 s, (f) t = 2 s, (g) t = 4 s, (11) t = 7 s.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated pa.~cles
359
of up to 4% between the solutions obtained enforcing P(t!4.1, 0)= 1 and those with no boundary condition enforced at the points (4-1, 0). Physically it can be argued in pro of either situation; however, from a numerical point of view, in the limit as ,.'ix, Ay ~ 0 the condition P(tl+-l, 0 ) = 1 should be enforced. The latter again exhibits an oscillatory behavior localized about the points (4.1, 0); however, as was the case before, this oscillation disappears very rapidly (within two nodes) and did not cause stability problems. The results presente6 in this work have been obtained in the Control Data Corporation C Y B E R 175 machines at the University of Illinois and the University of Arizona. A typical calculation in the 196 element mesh was obtained in 2.5 CPU s per time step, this time being by no means optimal.
8. Conclusions The numerical solution of degenerate elliptic equations has been recognized as a difficult task in computational methods, and the only reported results in the literature are those of Franklin and Rodemich [5] for the steady-state equations (5)-(6). The numerical solution of the transient equations (2)-(4) presented in this work is the first of its type, and in the absence of any analytical solutions to this problem provides a,l important tool for the analysis of a variety of problems, many arising in stochastic processes, which up to now have not been possible to solve. Examples of these can be found in [1. 18, 19]. An important application of the present algorithm lies in the solution for the probability of failure of the single degree of freedom linear oscillator [1. 2]. Solutions to this problem will provide valuable insight into more complex situations in structural reliability. Work is currently in progress in this area and the results will be reported in the near future.
Appendix A Consider a randomly accelerated particle whose position .\'(t) at lime t ,,atistics the stochastic differential equation =
n(t).
t>o.
with initial conditions x ( 0 ) = x, where
y.
n(t) is a stationary Gaussian white noise [1], with mean and covariance E[n(t)]=O,
E[n(t)n(t + r ) ] = 2,rSoS(r).
tA~)
So is the magnitude of the spectral density of the excitation, and S ( ) is the Dirac delta function. The system is linear, and its responses (X, ,Y) are the components of a M;~rkov vector. Thus, the theory of Markov processes applies.
360
L.A. Bergman, J.C. Heinrich. Galerkin and randomly accelerated particles
We start with the equation of continuity, the Chapman-Kolmogorov-Smoluchowski equation
f(~,~;tlx, y;to)= f:~f:®f(~,~;tl£,~;tl)f(~,~;t~lx, y;to)d£d~, to
at = -
a~
o~ ~'
(A.5)
or its adjoint, the backward Kolmogorov equation
--0% = y O-X + qrSo 0y 2
(A.6)
both subject to the initial condition
limf(£,9;tl x, y; to)--8(~- x)8(9- y),
(A.7)
t -~ to
and where af/at = -of/ato
(A.S)
.
When the system described by (A.1)-(A.3) is investigated for the first passage probability, a similar theory can be developed. Consider the 'safe' region. n - {(x, Y)[ Ix[ ~ 1, lyl < oo},
(A.9)
as shown in fig. 1. Eq. (A.5) now governs the evolution of the transition probability density function f(~, ~; t[(x, y; t0) N/~), where f d~ d~ is the probability that the particle be within states (£, ~) and (~ + d~, ~ + d~) at time t given that the initial state is (x, y) and that the particle has remained within t h e safe region f~ over the interval [t0, t]. The boundary conditions associated with (A.5) must now reflect the presence of absorbing boundaries at = -+1 in the phase plane. Defining the reliability function
F(tlx, y;to)=
f(£.~;tl(x,y;to)Nn)d~d~
(An0)
1
we see that we may substitute (A.8) into (A.6) and integrate term by term in accordance with
L.A. Bergman, Z C. Heinrich, Galerkin and randomly accelerated particles
361
(A.10) by interchanging the integration and differentiation operations. This yields an initialboundary value problem.
OF c~F OeF Ot = y - ~ + ¢rSo tgy2,
(A.11)
y<0,
F ( t ] - l , y; to)= 0, F(t] 1, y; to) = 0,
~y>0,
F(t l x, y ; to)-,O
as ]Yl--' oc,
F(tolX, y; to)= 1
for all (x, y).
(A.12) (A.13)
The probability density function of the first passage time distribution is given in terms of the reliability function as 0 g(tlx, y; to)= - ~ F(t l x, y; to).
(A.14)
The first moment of the first passage time distribution. T(x, y), can be expressed as (A.15)
T(x, y)= f,,] tg(t [ x, y; to)dt.
Multiplying each term of ( A . I I ) b y t and integrating term by term in accordance with (A.14) and (A.15) leads to the Pontriagin-Vitt equation
OT
~2T
(A.16)
- 1 = y - ~ + ¢rSo - ~ r .
T(-1, y ) = 0
y<0
t
T(1, y)= 0,
y > 0,
[
T(x, y)--, 0
as lyl--, oo.
(A.17)
Finally, recognizing that the probability of failure is the complement of the reliability function.
P(tlx, y; to) = 1 - F(tlx, y; to)
(A. 18)
an initial boundary value problem is easily derived for the probability of failure
OP OP ozP Ot = y O-x + ¢rSo 0y2 ,
(A.19)
eOl-l,y;to)=l,
y
1
P ( t l 1, y; to)- 1,
y >0,
t
P(tlx, y; to)-- 1 P(tolx, y;to)=O
as lyl -~ oo,
J
!
for all (x, y).
(A.20) (A.21)
362
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
Acknowledgement The authors are grateful for the support received from their respective departments for providing the necessary resources and computer time to perform this work. The assistance of Mr. Eric Austin in the preparation of graphical results of the analyses performed is also gratefully acknowledged.
References
[11 [21 [31 [41
Y.K. Lin, Probabilistic Theory of Structural Dynamics (McGraw-Hill, New York, 1967). J.N. Yang and M. Shinozouka, First-passage time problem, J. Acoust. Soc. Amer. 47 (1970) 393-394. S.H. Crandall, Fi~t crossing probabilities of the linear oscillator, J. Sound and Vibration 12 (1970) 285--299. L.A. Bergman and J.C. Heinrich, On the moments of time to first passage of the linear oscillator, Earthq. Engrg. Struct. Dyn. 9 (1981) 197-204. [51 J.N. Franklin and E.R. Rodemich, Numerical analysis of an elliptic-parabolic partial differential equation, SIAM J. Numer. Anal. 5 (1968) 680-716. [6] L.A. Bergman, Moments of time to first passage of the linear oscillator, Doctoral Thesis, Department of Civil Engineering, Case Western Reserve University, Cleveland, OH (1979). [7] J.R. Cannon and C.D. Hill, A finite-difference method for degenerate elliptic-parabolic equations, SIAM J. Numer. Anal. 5 (1968) 211-218. [81 R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, New York, 1967). [91 M. Abramowitz and I.A. Stegun, eds., Harldbook of Mathematical Functions (National Bureau of Standards, Washington, D.C., 1964). [1o; LC. Heinrich, P.S. Huyakorn, O.C. Zienkiewicz and A.R. Mitchell, An 'upwind' finite element scheme for two-dimensional convective-transport equation, Interaat. J. Numer. Meths. Engrg. 11 (1977) 131-143. [lll J.C. Heinrich and O.C. Zienkiewicz, Quadratic finite element schemes for two-dimensional convective-transport problems, Internat. J. Numer. Meths. Engrg. 11 (197'7) 1831-1844. [121 O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 1977) 3rd ed. [13] L.A. Bergman and J.C. Heinrich, Solution of the Pontriagin-Vitt equation for the moments of time to first passage of the randomly accelerated particle by the finite element method, Internat. J. Numer. Meths. Engrg. 14 (1980) 1408-1412. [141 J.C. Heinrich, On Bubnov-Galerkin, Petrov-Galerkin and Added Balancing Dissipation (Department of Aerospace and Mechanical Engineering, The University of Arizona, 1980). [15] P. Hood, Frontal solution program for unsymmetric matrices, Internat. J. Numer. Meths. Engrg. 10 (1976) 379-399. [16] D.F. Grifliths and J. Lorenz, An analysis of the Petrov-Galerkin finite element method, Comput. Meths. Appl. Mech. Engrg. 14 (1978) 39-69. [171 D.C. Dowson and A. Wragg, Maxunum entropy distribution having prescribed first and second moments, IEEE Trans. Information Theory (September 1973) 689-693. [181 P. 3amet, Numerical solution of the Dirichlet problem for eUiptic-parabolic equations, SIAM J. Numer. Anal. 6 (1969) 458--466. [191 B.J. Matkowsky and Z. Schuss, The exit problem for randomly perturbed dynamical systems, SIAM J. Appl. Math. 33 (1977) 365-382. [2o1 G.E. Myers, Numerically-induced oscillation and stability characteristics of finite-element solutions to twodimensional heat-conduction transients, University of Wisconsin Engineering Experiment Station. Report No. 43 (August 1977).
PETROV-.GALERFAN FINITE ELEMENT SOLUTION FOR THE FIRST PASSAGE PROBABILITY AND MOMENTS OF FIRST PASSAGE TIME OF THE RANDOMLY ACCELERATED FREE PARTICLE L.A. BERGMAN Department of Theoretical and Applied Mechanics. University of Illinois. Urbana. IL 61801. U.S.A.
J.C. HEINRICH Department of Aerospace and Mechanical Engineering. The University of Arizona. Tucson. A Z 85721. U.S.A.
Received I() August 19~() R e v i ~ d manuscript received 1~ January 1981
A numerical solution to an initial boundary value problem governing the probability of failure of a randomly accelerated free particle is obtained using a Petrov-Gaderkin tinite element method. This direct solution is the first successful one. ,'rod no o~hers ha~.e bccn reported ia the literature. A solution of the Pontriagin-Vit.t equation for the t|mc to first passage of the particle is obtained first: in this case an analytical solution is available and used to e~aluate the numerical algorithm. Extensions to the solution of other stochastic differential equations, in p~:rt|cular th('Jse governing the probability of failure of the linear oscillator, and applications to ,,tructur;d dsnam~cs are discussed
I. Introduction
In recent years a number of probability based codes have been formulated as a ba,,~is for safety analysis and design, providing a consistent framework for the swtematic as,~essment of the reliability of a structure. This evolution has occurred as a result of high risk technological development in structural systems. The determination of the probability of failure of such a system during its expected life is. in the most general sense, the reliability problem. When failure is defined as the first excursion of a system beyond some limit, whether displacement, strain or stress, due to a random excitation, the problem is characteristically referred to as a first passage problem. The simplest problem in structural reliability is represented by the linear single degree-offreedom oscillator subjected to stationary Gaussian white noise excitation, for which the displacement and velocity responses form the components of a two dimensional Markov vector process [1]. Therefore, the evolution of the transition probability density function is governed by a Fokker-Planck equation and, alternately, by its adjoint the backward Kolmogorov equation. Although the latter leads to a well-posed initial-boundary value problem for the system probability of failure allowing the use of direct numerical integration procedures, prior efforts 0045-7825/81/0000--0000/$02.50 © 1981 North-Holland
346
L.A. Bergman, 2.C. Heinrich, Galerkin and randomly accelerated particles
in this direction have been unsuccessful [2]. This has led to the formulation of a variety of approximations which are discussed in [1, 3, 4] and references therein. In the present work, we study a limiting case of the linear oscillator wherein the spring and dashpot are removed leaving a randomly accelerated free particle. The resulting initial boundary value problem contains all of the numerical difficulties of the more general problem with one fewer parameter. Furthermore, an analytical solution exists [5] for the PontriaginVitt equation governing the first moment of time to first passage of the free particle, permitting an assessment of the accuracy of the algorithm. A numerical solution of the transient equation for the free particle is obtained using a Petrov-Galerkin finite element method. This solution is, to the authors' knowledge, the first of its kind. The results presented here show that the method constitutes a powerful tool for the solution of a variety of problems arising in stochastic processes for which solutions have not before been available.
2. Problem description We consider a randomly accelerated particle which satisfies the stochastic differential equation
J~(t) = n(t) t > 0 ,
(1)
where X(t) is a random displacement process, n is a stationary Gaussian white excitation, the dot represents differentiation with respect to time, and with initial position and velocity X(0) = x, X(0) = y, over a 'safe' spatial domain O = {(x, y)/lxl- 1, lyl < o0}, which is shown in fig. 1. Let P(tlx, y) be the probability of failure of the particle to remain in the safe region, having started at (x, y), and without having previously failed and re-entered the safe domain. The function P then satisfies the initial boundary value problem c3P c3P c~2p Ot = y - ~ + ¢rSo 0y2 ,
(2)
with initial and boundary conditions
POI x, y) = 0
for (x, y) E n ,
P(t I-1, y)= 1
for y < 0 ,
P(t[ 1, y ) - 1 P(tl x, y ) ~ 1
for y > O,
(3)
(4)
as lyl--' oo,
where So is the magnitude of the spectral density of the excitation. A characteristic of this typ~ of equation is that no boundary condition is prescribed on part of the boundary, in. this c~,~ along the boundaries x = - 1 , y > 0 and x = 1, y < 0. A short account of the derivation of eqs
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
347
T=O p=l
. _.Y
T=O
p=t
(-I,O)
.,
(I,0)
(o, o)
x
sT p } FREE
T--O p=l
I I
i
-y T=O. p.'I
Fig. 1. Spatial domain and boundary conditions for cqs. (2) anti (¢i).
(2)-(4) is given in Appendix A. For further details the reader is referred to [1, 2.6]. it has been shown [2] that eqs. (2)-(4)constitute a well-posed problem: however to the author's knowledge no direct solution either analytical or numerical has been reported. Defining the function T(x, y) to be the expected value of the time to first passage of the randomly accelerated particle (i.e., the expected value for the first time the particle reaches a position IX(t)l = 1), T(x, y) is then the solution of the boundary value problem (see Appendix
A) 02T
OT
with boundary conditions
T(-l,y)=O
fory
T(I, y)= 0
for y > O,
Ttx, y)--, 0
as ly[--, 00.
(¢,)
Eq. (5) is known as the Pontriagin-Vitt equation for the first moment of time to first passage of the randomly accelerated particle. Problem (5)-(6) was solved analytically and numerically by Franklin and Roclemich [51 for
348
L.A. Bergman, J.C. Heinrich, Oalerkin and randomly acceleratedparticles
the case ¢rSo= 0.5. In their numerical calculations they used two finite difference methods previously developed in [7, 8]. However, their numerical results were not satisfactory, and they only obtained reasonable accuracy by averaging the nodal solutions in one of their methods. The analyst faces two basic difficulties when attempting the numerical solution of prob)em (5)-(6). The first stems from the character of the operator, Eqs. (2) and (5) fall into the general category of degenerate-elliptic or elliptic-parabolic operators which are positive semidefinite. In fact, eq. (5) has the form of a heat equation with time x and temperature T propagating from the right when y > 0 and from the left when y < 0. This requires that the numerical scheme be capable of taking the direction of the initial velocity y into account and producing discretizations that are biased accordingly. The oecond difficulty comes from the form of the solution of problem (5)-(6). Franklin and Rodemich [5] obtained an analytical solution that along y = 0 takes the form
T(x, 0 ) =
x2)~'6[F(-1/3,1, 7/6,½(1- x))+ F(-1/3, 1, 7/6, ½(1+ x))], C--3~"6/(2F(1/3)), F is the hypergeometric function and F the gamma C(1-
(7)
where function [9]. Clearly, the first partial derivatives of T with respect to x are unbounded at x = +_1, which causes additional computational difficulties and raises questions about convergence of the numerical algorithms. These will be discussed in a later section. In what follows we will assume ~rS0= 0.5, define a Petrov-Galerkin finite element method for the numerical solution of problem (5), (6), analyze its behavior and accuracy by comparison with the analytical solution (7) along the x-axis and discuss convergence. We then extend the algorithm to the transient equation (2), which is solved, to our knowledge, for the first time and evaluate the results.
3. Numerical solution of the Pontriagin-Vitt equation The partial differential equations (2) and (5) can be classified in general as convective diffusion equations in a highly anisotropic diffusive media. The difficulties related to the numerical solution of such equations when the first derivative (convective) terms are dominant led to the development of Petrov-Galerkin finite element methods by Heinrich et al. [10, 11] which are used in the context of the present problem. For reasons to be discussed in the next section, we restrict ourselves to a bounded rectangle ~ r = {(x, y)I Ixl-< 1, lYl-< Y} where Y is a large number (see fig. 1). The asymptotic boundary condition is prescribed at lYl = Y- Following a standard weighted residuals formulation [12], we replace the problem of finding the solution T(x, y) of the differential equation (5) over ~2~, by that of finding a function T(x, y) such that
fa Y{- ~IOWOT ~ ~yy 4 - y W ~OT} - d x d y - - fa
Wdxdy.
(8)
Y
T(x, y) is then
approximated by T(x, y)----Zk Nk(x, y)Tk, where Nk(x, y) are the well known bilinear quadrilateral basis functions and Tk the approximation to T at node k. For a particular element in isoparametric coordinates (~, ,/) the basis functions are given by
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
N'(~:, r/)= ¼(1-
349
~e)(1 - I,/)
N2(se, n)= I(1 + ~e)(1 - rl) N3(se, ~)= ¼(1 + ~e)(1+ n)
-1<_~, r/_< 1.
(9)
N'(~, n) = I(1 - O 0 + n) Following Heinrich et al. [10], the weighting functions are chosen to be
W2(~, 7/)= N2(g, 7/)(1 + 3 a ( 1 - ~:)),
W'(se, n) = N'(~:, "r/)(1 --~a(1 + ~)), W3(6 r/)= N3(se, r/)(1 +-~a(] - se)),
W4(~, 1~) -- N 4(~:, 7/)(1- ~a(1 + ~:)) ,
(10)
where t~ = 1 if y < 0, 0 if y = 0 and - 1 if y > 0. Fig. 2 shows the element notation together with two of the above weighting functions for the case t~ = - 1 , that is, for a node in the upper half plane. The other two are mirror images of the ones shown about the line ~, = 1. For t~ = 1 the weights are mirror images about the line ~ = - 1 . The effect of the special weighting functions in this numerical algorithm is better understood if we examine the difference equations generated for each node. We consider a rectangular mesh with nodes (xi, yj) at which T(x, yj)= T~. For simplicity we also approximate the coefficient y of the x-derivative by yj in the equation for node (x, y,). Under these conditions we obtain: (a) For y > 0 and a = - 1:
+
+ ~
+(-2Ax-2ydty2)T~+ - T + 2 y ~ y 2 +
y,A~22)T~-.~+\4 !
T~.,+ -
+Y/ 2 ] ~'~=-3Ax'~y2" N4
w4 4
Nil WI
T',:'~+Ax-y,
T',*' (II)
N3 w3
2NZ - g W2 '
I/*
.~
\ a=-I
350
L.A.
(b) For y
=
Bersman, ZC. Hein6ch, Galerkin and randomly accelerated particles
0 and a
=
O:
(12)
+ (~)T~+-I + (Ax)T~+' + ( ~ ) T { $ 1 = - 3 A x A y 2 •
(c) For y < 0 and a = 1 (here yj = [yj]):
* (--2Ax-- 2yjAy2)T~"['I-~)T~+I'Jt'(5~X3c--~J~2~erJ+'21i-1"~" ( Ax-y]~--~-2 )T,'+I + (-- T AX~'r'/+' ! " ÷ ' = -3AxAy 2
(13)
In fig. 3a we further illustrate the nine node stencil for the case y > 0 and a = - 1. It can be seen that the Petrov-Galerkin finite element method not only approximates the first derivative term in a backward way with respect to the velocity direction bu~ it also shifts the weighting of the second order term in the x-direction in the nine node stencil, still maintaining the symmetry of the approximation in the y-direction. This property of Petrov--Galerkin schemes, to date not well understood, appears to be a key factor in the algorithm's success. Furthermore, in the presence of
Ax-yj AY22
O
O
(Xt_l,Yj+ 1)
(xi, Y:]+I)
(x.I+I,YJ+ 1)
~x --
-2~x-2y ~y2 Oj
_ 5~x4+ 2yjAy2
- 7~x4-2YJAy2
O (xi+I, Yj)
(x i ,
(Xt_l,Yj) ~x ---~ O (xi_ 1, Yj_l)
(xt, Yj) 4y2 ~x-yj 2 O (xi, Yj_I)
a
5Ax8+ YJAy22
7Ax _Ay2 s "Yl 2
- -'8Az
O
8
YJ 2 O
(xi+1,
Yj_I)
5Ax8+ YJAy22
O'
0
(xl] YJ+I)
(xi+1,
/j+l)
_ 5~x4j 2YJAyz
O
O (xi+ 1, Yj)
Yj)
7Ax Ay2 8 - YJ 2 O (xi, Yj_l)
5Ax [ 8
~y2
YJ 2
O
(xi+I, Yj_I) b
Fig. 3. Differences stencils for Petrov-Galerkin Finite Element Method; (a) full molecule for interior node, (b) molecule for a node at the free left-hand side boundary.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
351
non-constant source terms, as is the case if higher moments are to be computed using the Pontriagin-Vitt equation [13, 4], it can be easily proved that the limiting one-dimensional case dT/dx = Q,
(14)
T(O)= 0
is exactly integrated when O is linear [14]. Any other existing backward differencing scheme applied to equation (14) will give a numerically diffused solution in the presence of a linear source term O, the difference in the Petrov-Galerkin scheme being that O is weighted consistently, with the consequent re-distribution of the nodal contributions of the source term to the right hand side vector. In fig. 3b we show the contributions to the stiffness matrix for the case of a node in the part of the boundary at x = - 1 where no boundary condition is imposed. It is clear, since the derivative of T with respect to x is approximated by a fully backward difference with respect to the velocity, that the numerical scheme does not impose any kind of boundary condition there. In fact, the present scheme is free of any kind of natural boundary condition. The resulting discretization is assembled in the form (15)
KT= B,
where
k° =
v
1 8W' aN__~ j + yW~
2 8y
c~y
Ox J
dx dy
and
b,
~,
dx d y .
This system of linear equations is solved by a direct frontal elimination procedure for unsymmetric matrices [15].
4. Convergence of the Petrov--Galerkin method Several questions arise concerning convergence of the numerical algorithm presented in the previous section. At this point, it is not within our capabilities to present a proof for the convergence of the method. In fact, none of the, existing theorems appears to be applicable in this case. Our arguments will be later substantiated by the numerical evidence. The fact that no second derivative term with respect to x appears in eq. (5) implies that c~ does not go to zero as the mesh size decreases. This makes the results of Griffiths and Lorenz [16] not applicable to our equations. A straightforward truncation error analysis applied to eqs. (11)-(13) indicates that the truncation error is O ( A x + dy:) as expected. However the leading terms depend on the higher order derivatives ~3Y/o~x~y 2 ,
02T/~x 2
and
8"T/Sy 4 .
Examination of equation (7) reveals that along the x-axis the solution T is the product ~,f an even function of x by ( l - x 2 ) TM, and the first partial derivatives with respect to x are unbounded in a neighborhood of x = _+1, y = 0. Furthermore the square of the first derivative
352
L.A. Bergman, .I.C. Heinrich, Galerkin and randomly accelerated particles
is non-integrable. As a consequence, in this particular case the truncation error in (11), (13) is unbounded and from the Petrov-Galerkin point of view we are violating the basic requirement of the solution lying in the space of functions with a finite energy. This was resolved by Franklin and Rodemich [5] by constructing a family of similar functions T(x, y, a) such that (a) T(x, y, a) is smooth for Ix[--< 1, [y[-< Y, if a > 1, (b) T(x, y, a)--, T(x, y) as a -- 1, and (c) LT(x, y, a ) ~ LT(x, y) as a -~ 1, where L =-y(O/Ox)-~(0Z/0y:). These functions essentially shift the singularity to the points x = +-a, y = 0 with ial > 1 as a consequence though no convergence rates will be obtainable. For more details see Franklin and Rodemich [5]. The fact that the domain is unbounded in the y-direction presents no difficulty from the finite element point of view since asymptotically the solution has the behavior 1-x Y
T(x, y)
as y-,0o,
T(x, y ) ~ - l~+ x Y
as y - , - 0 o ,
(16)
and this is square integrable over the infinite region. However, in order to avoid dealing with an infinite region and unnecessary complications introduced by the use of 'infinite elements' (see e.g. [12]) we choose to truncate the region at a far enough distance Y and accept introducing yet another small error in the manner of Franklin and Rodemich [5]. It will be clear later from the numerical results that the error introduced by this approximation is extremely localized and has no significant effect even if homogeneous boundary conditions are imposed at the truncated boundary for fairly small values of Y.
5. Numerical solution of the transient equation
We now turn to the solution of problem (2)-(4). The spatial operator remains unmodified from the previous case so now the Petrov-Galerkin finite element method is applied to eq. (5) with the right hand side replaced by the time derivative. This yields the semi-discretized system M//-KP
= 0,
(17)
where the stiffness matrix K is the same as in eq. (15), P is the vector of nodal approximations Pk to P(tlxE, Yk) and/~ is its first time derivative. The mass matrix M is now given by
m° = fn WiNS dx dy.
(18)
Y
Following the notation of section 3, the resulting assembled mass matrix for a typical node (xi, Ys) is shown in fig. 4a. We point out that to make it compatible with the expressions for the stiffness matrix of fig. 3 the terms in the mass matrix must be multiplied by 3Ay. An interesting situation arises when a boundary node with no boundary condition applied is
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles AxAyr~ a 12 (3 + ~)
AxAy.I
AxAy
O'
Yj+l)
9 O
~C i 12 "3 + 0
a 2)
4AxAy
I
AxAy.2
-12 ( 3 '0
0
(xi_I, Yj_I)
(x t , Yj_I)
I,
AxAy:! a O
0
(xi+I, Yj)
AxAy.1
(Xi+l~ Yj+I)
3 t~ - ~)
o
AxAv 9
O
[ Y]+I)
(xi'
(xi[ YJ)
A:'~y.l a
ct
O
O
(xi+I, Yj+I)
9 O
(xi_1, Yj)
AxAy.2
(xi, Yj+!)
(xi-Ii
nxAy.t a 3 ~3 + ~) O
a
n 2)
(xi+1, Yj_!)
a
(xi+I, Yj)
(xii YJ)
AxAy.2 .,~ 12 (3 - 2)
AxAy I - ~) 12 (3
0 (x.i ,
353
0
Yj_!)
(xi+I, Yj_I) b
Fig. 4. Mass matrix for Petrov-Galerkin Finite Element Method; (a) interior node. (b) node at the free left-hand side boundary.
reached, as shown in fig. 4b for the case y > 0. Since a = - 1 in this case, the method will increase the mass at these boundaries. This should result in a stabilizing effect there, overall conservation being achieved by subtraction of the same amount at the opposite boundary. This, however, indicates that a mesh symmetric about the y-axis should be used. It is clear in fig. 4 that the Petrov-Galerkin method redistributes the masses of the non-central nodes to the point of assigning negative masses to those downstream with respect to the velocity direction, while keeping the symmetry in the y-direction. Eq. (17)is further discretized using the theta method [12]: that is [M-
At(1 - 0)Kle
"+' = [ M -
AtOK]P'.
(l.)
where 0 <-0-< 1, n denotes the current time and n + 1 the time to advance to. For 0 = 0.5 the second order in time Crank-Nicolson marching scheme is obtained, and for 0 = 0 the fully implicit forward marching method, both unconditionally stable. These two methods have been used for the time discretization with time steps varying in the range 0 . 0 ~, ~ A t ~ O.2 s. T h e critical time step, beyond which a fully explicit (0 = 1.0) .scheme would be unstabie, was estimated in accordance with Myers [20] to be 0.025 s. An upper bound for the time steps prior to excessive loss of accuracy or the onset of oscillations was found to be dependent on the problem and in particular on the mesh.
6. Numerical results for the Pontriagin-Vitt equation Numerical solutions to problem (5), (6) have been obtained fol a var:~ety of situations to assess the effect of the mesh, boundary conditions and the artificial boundary y = ± Y.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
354
The results indicate that the boundary conditions at y = - Y have a very localized effect, since we are interested mainly in the results about y = 0 where the maximum gradients occur. We find that for a value of Y as small as 5 no difference can be detected in the solution along y = 0 when the asymptotic boundary conditions (16) are replaced by homogeneous conditions. A second question is how the imposition of a homogeneous boundary condition at the points (+-1, 0) affects the solution (Franklin and Rodemich [5] did not impose T(+I, 0)= 0 in one of their finite differences schemes). It was found that although the imposition of T(+__l, 0)= 0 causes some oscillation close to those nodes, this oscillation disappears rapidly (within two nodes) and can be reduced considerably by adequate mesh refinement in the neighborhood of the singular points. The relaxation of this condition on the other hand results in a small overall loss of accuracy, which is minimized through mesh refinement, however no oscillations are observed in this case. From the above it also becomes clear that the best choice of mesh is an irregular one refined in the neighborhood of the points (x, y ) = ( ± l , 0), regardless of the boundary condition strategy. Here we present the results obtained with a 20 by 20 element mesh, with nodes symmetrically distributed about the x and y axis with coordinates at x = 0.0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.975, 1.0 and at y = 0.0, 0.2, 0.4, 0.7, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0. The artificial boundary is at Y = 5, and homogeneous boundary conditions were imposed there. In table 1 we show results along the x-axis obtained with and without imposing the conditions T(+_l, 0)= 0. These are compared with the analytic solution and with the finite difference solutions of Franklin and Rodemich [5] using an 11 by 21 nodes mesh with Y = 20. Fig. 5 shows the finite element solution obtained imposing the condition T(_+I, 0)= 0, for initial velocities y = 0.0, 0.4, 1.0 and 2.0; and in fig. 6 we show the three dimensional surface for the entire solution. The agreement of the finite element calculations with the analytical solution is excellent (exact to three significant figures) as shown in table 1. We point out again that refinement about the points (x, y) = (±1, 0) plays a critical role in the accuracy of the solution, and nothing can be said about convergence rates there. Table 1 Analytical and numerical results for the Pontriagin-Vitt equation X
EXACT
FEI a
FE2 b
FD1 c
FD2 ~
X
EXACT
FD3 "
0.0 0.2 0.4 0.6 0.8 0.9 1.0
2.312 2.289 2.216 2.076 1.813 1.580 0.000
2.314 2.290 2.217 2.077 1.826 1.578 0.000
2.317 2.29g 2.220 2.080 1.828 1.598 0.762
2.151 2.125 2.042 1.878 1.554 -0.000
2.475 2.115 2.409 1.850 2.171 m
0.000 0.182 0,364 0.545 0.727 0.909
2.312 2.293 2.233 2.123 1.933 1.550
2.319 2.268 2.246 2.088 1.980 1.432
"FE1 = Finite elements with T ( - 1 , 0) = 0. Finite elements with T(4-1, 0) free. CFD! = Biased finite differences. dFD2 = Biflanced finite differences (no averaging). t "FD3 = Balanced finite differences (averaged). bFE2 =
Ref. [5].
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
355
~ " , , ~ y =0.0 2.0 I..-
O0 INITIAL DISPLACEMENT x
Fig. 5. Expected values of time to first passage T vs. initial displacements for initial velocities y = O.O, 0.4. 1.0 and 2.0.
x
o
•
v
,
•
~
Fig. 6. Three dimensional surface for the expected value of time to first passage T
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
356
7. N u m e r i c a l r e s u l t s f o r the transient e q u a t i o n
Considering the previous experience in the solution of the Pontriagin-Vitt equation (5), numerical experiments for the solution of the initial boundary value problem (2) were restricted to the case Y = 5. We present results of computations performed in a 196 element irregular mesh with nodes symmetrically distributed about the x and y axis and coordinates Ix[ = 0.0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1.0 and [y[ = 0.0, 1.0, 2.0, 3.0, 3.5, 4.0, 4.5, 5.0. Table 2 shows results obtained using both the Crank-Nicolson and fully implicit methods, for the cases P(t[_+1,0)= 1 and P(ti -_.1, 0) free. In the same table we show some values for T(0, 0) obtained by calculation of the first moment of the numerical solution i.e.
Ktg(tlx, y;n)dt,
T(O,O)=
(20)
w h e r e g ( t l x , y ; D ) is t h e p r o b a b i l i t y density of t i m e to first passage (see A p p e n d i x A). C o m p a r i s o n with t h e exact ~Jiation s h o w s less t h a n 5 % error. Fig, 7 s h o w s t h e probability of failure at t h e p o i n t (0, 0) versus time, a n d it is c o m p a r e d with t h e results o b t a i n e d using a t w o m o m e n t m a x i m u m e n t r o p y distribution [17] with o r d i n a r y m o m e n t s nl = 2.295 a n d n2 = 7,29, a n d an e x p o n e n t i a l (one p a r a m e t e r ) distribution [6, 13]. A n excellent overall a g r e e m e n t can b e
Table 2 Evolution of the probability of failure at the origin P(z [0, 0); Petrow-Galerkin finite element solution
P(t tO, O) t (s) 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
1"
2b
3c
0 . 0 0 0 0.000 0 . 0 8 9 0.062 0.462 0.484 0.707 0.722 0 . 8 3 7 0.848 0 . 9 1 0 0.917 0 , 9 5 0 0.955 0 . 9 7 3 0.97.6 0 . 9 8 5 0.987 0.992 0.993 0 . 9 9 5 0.996
T(0, 0) [Exact = 2.312]
"At = 0.2, 6' -- 0.0 bat = 0.1, 0 = 0.0 ~At = 0.05. 0 = 0.0
0.000 0.042 0.497 0.729 0.853 0.921 0.958 0.977 0.988 0.993 0.996
A.d
0.000 0.000 0.495 0.728 0.85J, 0,9"23 0.959 0.978 0.989 0.994 0.997.
2.483
P(t [ +--1,0) Free.
dAt = 0-1, 0 = 0-5 J
~At = 0.1, O = 0.0 / fat = O.1, 0 = 0.5 f
P(t l +-1,O)= 1.
5e
0.000 0.065 0.506 0.748 0.870 0,933 0.966 0.983 0.991 0.995 0.998
6f
0.000 0.000 0.517 0.756 0.877 0.939 0.969 0.985 0.992 0.996 0.998 2.405
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles m
1.0
..-- -"
f
//
~
""
" -
. . . . . . . . . . . . . . . .
357 , ° .
o°'°"
.....-""
°°°°°
0.8
...~r/
n~ _J
-
0.6
..""/y •"
h
0 )-
_J ,~ .ff
ID 0Iz
..."
0.4
-
/
77
.""
FEM
0.2
O0 0
2
4
6
8
~0
TIME (SEC)
Fig. 7. Probability of failure P vs. time for initial conditions of zero displacement and velocity.
IO
- -
08 ;
iAL DISPLACEMENT AND vELOCF~;Y
06
(~,y} = ( - 0 6 , 0
.. " 04
-
:(-04,
00,~
/
"
I
;
/
.
I
."
(=,y}
/
.
0,2
.
O)
/ //
oo
o
2
4
6
8
tO
T I M E (SEC)
Fig. 8. Probability of failure P vs. time for three different nonhomogeneous initial conditions of displacement and velocity.
358
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
observed in fig. 7 with a significant correction to the approximations shown by the finite elemev.t solution for lower times 0 < t < 1.5. In fig. 8 similar curves are shown for the points (0, 2), (-0.4, 0), (-0.6, 0). These plots come from run Number 3 in table 2. Fig. 9 shows the evolution of the three dimensional surface with time. After seven seconds the change in the value P(t 10, 0) become less than 0.2% every 0.1 second, and changes are no longer observable. In general few difficulties were encountered in the solution of the present problem. The initial time steps were subject to some oscillatory behavior due to the singularity in the boundary conditions. The fully implicit method (0 = 0) shows less disturbances than the Crank-Nicolson method (0 = 0.5) as would be expected, although the latter yields what appears to be a more accurate solution at later times. Also, in table 2 we observe differences
o o,>c, ~''~'~
a
t
=0 l se¢
t = 0 0 sec
c
t :05sec
t --I O s e c
t= t=2Osec.
f = 4.0see
t ='/'.0 $eC.
Fig. 9. T i m e e v o l u t i o n of t h e t h r e e d i m e n s i o n a l p r o b a b i l i t y of f a i l u r e surface; (a) initial c o n d i t i o n s , (b) t = 0.1 s, (c) t = 0.5 s, (d) t = 1.0 s, (e) t = 1.5 s, (f) t = 2 s, (g) t = 4 s, (11) t = 7 s.
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated pa.~cles
359
of up to 4% between the solutions obtained enforcing P(t!4.1, 0)= 1 and those with no boundary condition enforced at the points (4-1, 0). Physically it can be argued in pro of either situation; however, from a numerical point of view, in the limit as ,.'ix, Ay ~ 0 the condition P(tl+-l, 0 ) = 1 should be enforced. The latter again exhibits an oscillatory behavior localized about the points (4.1, 0); however, as was the case before, this oscillation disappears very rapidly (within two nodes) and did not cause stability problems. The results presente6 in this work have been obtained in the Control Data Corporation C Y B E R 175 machines at the University of Illinois and the University of Arizona. A typical calculation in the 196 element mesh was obtained in 2.5 CPU s per time step, this time being by no means optimal.
8. Conclusions The numerical solution of degenerate elliptic equations has been recognized as a difficult task in computational methods, and the only reported results in the literature are those of Franklin and Rodemich [5] for the steady-state equations (5)-(6). The numerical solution of the transient equations (2)-(4) presented in this work is the first of its type, and in the absence of any analytical solutions to this problem provides a,l important tool for the analysis of a variety of problems, many arising in stochastic processes, which up to now have not been possible to solve. Examples of these can be found in [1. 18, 19]. An important application of the present algorithm lies in the solution for the probability of failure of the single degree of freedom linear oscillator [1. 2]. Solutions to this problem will provide valuable insight into more complex situations in structural reliability. Work is currently in progress in this area and the results will be reported in the near future.
Appendix A Consider a randomly accelerated particle whose position .\'(t) at lime t ,,atistics the stochastic differential equation =
n(t).
t>o.
with initial conditions x ( 0 ) = x, where
y.
n(t) is a stationary Gaussian white noise [1], with mean and covariance E[n(t)]=O,
E[n(t)n(t + r ) ] = 2,rSoS(r).
tA~)
So is the magnitude of the spectral density of the excitation, and S ( ) is the Dirac delta function. The system is linear, and its responses (X, ,Y) are the components of a M;~rkov vector. Thus, the theory of Markov processes applies.
360
L.A. Bergman, J.C. Heinrich. Galerkin and randomly accelerated particles
We start with the equation of continuity, the Chapman-Kolmogorov-Smoluchowski equation
f(~,~;tlx, y;to)= f:~f:®f(~,~;tl£,~;tl)f(~,~;t~lx, y;to)d£d~, to
at = -
a~
o~ ~'
(A.5)
or its adjoint, the backward Kolmogorov equation
--0% = y O-X + qrSo 0y 2
(A.6)
both subject to the initial condition
limf(£,9;tl x, y; to)--8(~- x)8(9- y),
(A.7)
t -~ to
and where af/at = -of/ato
(A.S)
.
When the system described by (A.1)-(A.3) is investigated for the first passage probability, a similar theory can be developed. Consider the 'safe' region. n - {(x, Y)[ Ix[ ~ 1, lyl < oo},
(A.9)
as shown in fig. 1. Eq. (A.5) now governs the evolution of the transition probability density function f(~, ~; t[(x, y; t0) N/~), where f d~ d~ is the probability that the particle be within states (£, ~) and (~ + d~, ~ + d~) at time t given that the initial state is (x, y) and that the particle has remained within t h e safe region f~ over the interval [t0, t]. The boundary conditions associated with (A.5) must now reflect the presence of absorbing boundaries at = -+1 in the phase plane. Defining the reliability function
F(tlx, y;to)=
f(£.~;tl(x,y;to)Nn)d~d~
(An0)
1
we see that we may substitute (A.8) into (A.6) and integrate term by term in accordance with
L.A. Bergman, Z C. Heinrich, Galerkin and randomly accelerated particles
361
(A.10) by interchanging the integration and differentiation operations. This yields an initialboundary value problem.
OF c~F OeF Ot = y - ~ + ¢rSo tgy2,
(A.11)
y<0,
F ( t ] - l , y; to)= 0, F(t] 1, y; to) = 0,
~y>0,
F(t l x, y ; to)-,O
as ]Yl--' oc,
F(tolX, y; to)= 1
for all (x, y).
(A.12) (A.13)
The probability density function of the first passage time distribution is given in terms of the reliability function as 0 g(tlx, y; to)= - ~ F(t l x, y; to).
(A.14)
The first moment of the first passage time distribution. T(x, y), can be expressed as (A.15)
T(x, y)= f,,] tg(t [ x, y; to)dt.
Multiplying each term of ( A . I I ) b y t and integrating term by term in accordance with (A.14) and (A.15) leads to the Pontriagin-Vitt equation
OT
~2T
(A.16)
- 1 = y - ~ + ¢rSo - ~ r .
T(-1, y ) = 0
y<0
t
T(1, y)= 0,
y > 0,
[
T(x, y)--, 0
as lyl--, oo.
(A.17)
Finally, recognizing that the probability of failure is the complement of the reliability function.
P(tlx, y; to) = 1 - F(tlx, y; to)
(A. 18)
an initial boundary value problem is easily derived for the probability of failure
OP OP ozP Ot = y O-x + ¢rSo 0y2 ,
(A.19)
eOl-l,y;to)=l,
y
1
P ( t l 1, y; to)- 1,
y >0,
t
P(tlx, y; to)-- 1 P(tolx, y;to)=O
as lyl -~ oo,
J
!
for all (x, y).
(A.20) (A.21)
362
L.A. Bergman, J.C. Heinrich, Galerkin and randomly accelerated particles
Acknowledgement The authors are grateful for the support received from their respective departments for providing the necessary resources and computer time to perform this work. The assistance of Mr. Eric Austin in the preparation of graphical results of the analyses performed is also gratefully acknowledged.
References
[11 [21 [31 [41
Y.K. Lin, Probabilistic Theory of Structural Dynamics (McGraw-Hill, New York, 1967). J.N. Yang and M. Shinozouka, First-passage time problem, J. Acoust. Soc. Amer. 47 (1970) 393-394. S.H. Crandall, Fi~t crossing probabilities of the linear oscillator, J. Sound and Vibration 12 (1970) 285--299. L.A. Bergman and J.C. Heinrich, On the moments of time to first passage of the linear oscillator, Earthq. Engrg. Struct. Dyn. 9 (1981) 197-204. [51 J.N. Franklin and E.R. Rodemich, Numerical analysis of an elliptic-parabolic partial differential equation, SIAM J. Numer. Anal. 5 (1968) 680-716. [6] L.A. Bergman, Moments of time to first passage of the linear oscillator, Doctoral Thesis, Department of Civil Engineering, Case Western Reserve University, Cleveland, OH (1979). [7] J.R. Cannon and C.D. Hill, A finite-difference method for degenerate elliptic-parabolic equations, SIAM J. Numer. Anal. 5 (1968) 211-218. [81 R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, New York, 1967). [91 M. Abramowitz and I.A. Stegun, eds., Harldbook of Mathematical Functions (National Bureau of Standards, Washington, D.C., 1964). [1o; LC. Heinrich, P.S. Huyakorn, O.C. Zienkiewicz and A.R. Mitchell, An 'upwind' finite element scheme for two-dimensional convective-transport equation, Interaat. J. Numer. Meths. Engrg. 11 (1977) 131-143. [lll J.C. Heinrich and O.C. Zienkiewicz, Quadratic finite element schemes for two-dimensional convective-transport problems, Internat. J. Numer. Meths. Engrg. 11 (197'7) 1831-1844. [121 O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, London, 1977) 3rd ed. [13] L.A. Bergman and J.C. Heinrich, Solution of the Pontriagin-Vitt equation for the moments of time to first passage of the randomly accelerated particle by the finite element method, Internat. J. Numer. Meths. Engrg. 14 (1980) 1408-1412. [141 J.C. Heinrich, On Bubnov-Galerkin, Petrov-Galerkin and Added Balancing Dissipation (Department of Aerospace and Mechanical Engineering, The University of Arizona, 1980). [15] P. Hood, Frontal solution program for unsymmetric matrices, Internat. J. Numer. Meths. Engrg. 10 (1976) 379-399. [16] D.F. Grifliths and J. Lorenz, An analysis of the Petrov-Galerkin finite element method, Comput. Meths. Appl. Mech. Engrg. 14 (1978) 39-69. [171 D.C. Dowson and A. Wragg, Maxunum entropy distribution having prescribed first and second moments, IEEE Trans. Information Theory (September 1973) 689-693. [181 P. 3amet, Numerical solution of the Dirichlet problem for eUiptic-parabolic equations, SIAM J. Numer. Anal. 6 (1969) 458--466. [191 B.J. Matkowsky and Z. Schuss, The exit problem for randomly perturbed dynamical systems, SIAM J. Appl. Math. 33 (1977) 365-382. [2o1 G.E. Myers, Numerically-induced oscillation and stability characteristics of finite-element solutions to twodimensional heat-conduction transients, University of Wisconsin Engineering Experiment Station. Report No. 43 (August 1977).