unnuun-nv- '93;-

0‘"...nun-.nuuauu

 

.,.;.........ou

 

 

 

 

l—
IES

’ liiilli\iillililg\\l\\i

 

LIBRARY
Michigan State

L University }

 

 

This is to certify that the
dissertation entitled
Some Numerical Methods For Singular
Diffusions Arising In Genetics

presented by

Nacer Eddine Abrouk

has been accepted towards fulfillment
of the requirements for

PhoDo degree in StatiStiCS

 

 

W4. 4. 5.“

Major professor

Date November 20, 1992

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

PLACE IN nETunN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

w . u
U 9 13%,
if: g. 11

2 8 12:94

, ‘1 ')
[X'gT
I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU In An Affirmative Action/Equal Opportunity Institution
cMMnmS-DJ

SOME NUMERICAL METHODS FOR SINGULAR
DIFFUSIONS ARISING IN GENETICS

By

Nacer Eddine Abrouk

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Statistics and Probability

1992

ABSTRACT

SOME NUMERICAL METHODS FOR SINGULAR
DIFFUSIONS ARISING IN GENETICS

By

Nacer Eddine Abrouk

In this dissertation, several different methods for approximating the laws of some
diffusions which are sOlutions of stochastic differential equations are considered.
Using Ito’s formula we derive an infinite system of differential equations for the
moments of the process. This system of ordinary differential equations is very
difficult to solve in most cases. A power series method is proposed for which
theoretical and numerical studies are made. Using the power series method, we
are able to approximate the moments of the diffusion under consideration and
to investigate its steady state behavior. The power series method leads to new
numerical algorithms for the study of the steady state analysis of the difl'usions
under consideration. The Gauss Galerkin method which is an existing method
originally introduced by D. A. Dawson in 1980 is also used. For the Gauss Galerkin
method we derive an alternative computational algorithm which has proved to be
efficient in practice. A detailed comparison between the power series method and
the Gauss Galerkin method is presented. These methods are then applied to some
stochastic models arising in population genetics. Numerical results relating to these
models are presented and discussed. Suggestions for future work and some open
problems, such as the approximation of the transition probability function, are

mentioned.

TO INGRID AND MICHAEL

ACKNOWLEDGMENTS

I wish to express my sincere thanks to Professors David Yen and Habib Salehi for
the suggestion of the problem and all the helpful direction in solving it. Their ad-
vice, their encouragement and the amount of time spent on discussions are greatly
appreciated. I would like to thank Professors Roy Erickson and Raoul LePage for
their time and interest. A special thanks goes to my wife, Ingrid, for her unwaver-

ing support and infinite patience.
Finally, I would like to thank the Department of Statistics and Probability and the

Department of Mathematics, at Michigan State University, for supporting me and

providing access to their computing laboratories.

iii

Contents

1 Introduction
1.1 Motivation ................................

1.2 An overview of the dissertation ....................

2 Statement of the Problem
2.1 Background Material from Probability ................
2.2 Background Material from Numerical Analysis ............
2.3 Statement of the Problem .......................

3 Numerical Methods
3.1 Assumptions and Derivation of the Moments System of Equations .
3.2 Summary of Previous Methods ....................
3.2.1 Moment Neglect Method ....................
3.2.2 Fixed Nodes Method ......................
3.2.3 Gauss-Galerkin Method ....................
3.3 The Power Series Method .......................

4 Some Models From Genetics
4.1 Description of the Various Models ...................
4.2 Solution of the System of Moments ..................
4.2.1 Moment Neglect Method ....................
4.2.2 Gauss-Galerkin Method ....................
4.2.3 Power Series Method ......................

iv

10

15
15
20
2O
22
24
29

38
38
40
4O
42
42

V

4.2.4 Comparison of the Power Series Method and the Gauss-Galerkin

Method ............................. 44
4.3 Steady State Analysis .......................... 45
4.4 Discussion ................................ 50
Appendix 53

Bibliography 71

Chapter 1

Introduction

1.1 Motivation

Let X = (X.;t 2 0) be a diffusion process arising from the stochastic differen-
tial equation: dX, = a(X,)dt + b(X,)dB, where X0 is given , B, is the standard
Brownian motion and the functions a and b are the drift and diffusion coefficients,
respectively. Such equations model many problems in biological, physical and engi-
neering sciences (for more details see Karlin and Taylor, Chapter 15). For a rigorous
mathematical treatment of theses models, we refer to Feller, where a review of the
theoretical work on diffusion processes as well as their boundary behavior is thor-
oughly carried out. Equations of this type play an important role in population
genetics. Wright (1931 and later publications) and Kimura derived the appropri-
ate discrete models in population genetics. The continuous counterparts of these
discrete models, which have been used extensively, are given by the stochastic dif-

ferential equation above (see Ethier and Kurtz).

Since, in general, the probability law of X“ t Z 0, is difficult to determine, an
important task is the problem of numerically approximating the probability law of
X;, for t Z 0. Because most population genetic models have singularities at the
boundary of the state space of the process X, it is even more difficult to approximate

the probability density (when it exists) of X,, for t 2 0, through the Fokker-Planck

2

equation for general drift and diffusion coefficients a and b. However, when the
coefficients a and b2 are polynomials, a system of equations for the moments of the

process X can be derived through the Ito formula (see Chapter 3 system 3.4).

If a is at most linear and b2 is at most quadratic, then this system is closed in the
sense that an explicit solution for the process moments is possible. In any other
case, this system of equations forms an ”open hierarchy” in which the equation
for the n-th moment may have terms involving the (n+1)-th moment or higher.
various truncation schemes have been proposed. The simplest and most frequently
used is the moment neglect method in which the (n+l)-th moment and higher mo-
ments are set equal to zero. In this way the system of moment equations reduces
to a closed system of ordinary differential equations. Unfortunately, as we shall
indicate below (see Chapter 3, section 3.2.1) the use of the moment neglect (or any
of its variants) can lead to misleading results. The Gauss-Galerkin method was
originally proposed by Dawson (1981) and then refined and improved by HajJar
far, and HajJafar, Salehi and Yen. It combines elements of the method of Gauss
quadratures and Galerkin approximation and yields an approximation in the sense
of weak convergence for the law of X t by discrete measures. It transforms the prob-
lem of solving the Fokker-Planck equation for the density of the law of X, (when it
exists) to one for the ”nodes” and ”weights” of the discrete measures (see Chapter
3, section 3.2.3). This method, when applicable, works for a more general class
of models which includes models where a and b2 are polynomials. In the latter
case, system 3.15 describing the dynamic behavior of the ”nodes” and ”weights”
can be viewed as a finite linear system of equations for the moments of the discrete
measure. The latter system is a truncation of the infinite linear system of ordinary
differential equations satisfied by the true moments of the process X (see Chapter
3, system 3.9). However, the linear system 3.9, when viewed as a system for the
moments of the discrete measure, is not closed in the sense that the equation for
the (2n-1)-th moment may involve the 2n-th moment or possibly higher moments.
Lemma 3.1 of Chapter 3 provides a procedure, called the Hankel truncation, to

close this system. The resulting system (system 3.17) is a nonlinear system of

3

2n equations for the moments of the discrete measure. When this system (system
3.17) has a solution on some t interval of the form [0,T] and T is a positive constant
independent of n, the Gauss-Galerkin method has proved to provide good approx-
imations in the sense of weak convergence of measures (see HajJ afar, Salehi and
Yen, 1992). Unfortunately, there are models in population genetics where T is a
function of n, causing the solution of system 3.17 to exist only locally with respect to
the variable t. However, numerically the solution seems to be global. Moreover, in

the Gauss-Galerkin method the initial distribution is required to be nondegenerate.

Our work in this dissertation consists mainly of the following contributions:

1. An improved numerical algorithm for the Gauss-Galerkin method when the

the coefficients a and b2 are polynomials (see Chapter 3).

2. A new approximation method, the power series method. The power series
method is an alternative to the Gauss-Galerkin method. This method yields a
way of replacing the infinite system for the moments of the process by a finite
system which is obtained by truncating the power series expansions of the
first few moments of the process (see Chapter 3, system 3.22). The resulting
system is a finite linear system which has a unique global solution and the
initial distribution is allowed to be degenerate. The power series method
yields a uniform (in the t variable) approximation for any j-th moment of
the process. However, the question of whether the power series method gives
rise to convergence in the sense of measures, as was the case for the Gauss-

Galerkin method, is not resolved.
3. An analysis of the steady state solution through the power series method.

4. A numerical comparison between some of the existing methods and the power

series method.

4

1.2 An overview of the dissertation

This dissertation is organized as follows. Chapter 2 contains background mate-
rial from both probability theory and numerical analysis, and the statement of the
problem. We begin with a review of the notions of diffusion processes, transition
probability functions and the associated semigroup. We then present some back-
ground materials from numerical analysis. This will be followed by a statement of

the problem.

We begin Chapter 3 with a detailed discussion of the assumptions needed. We
also show how the boundary conditions are important for the proper formulation
of the problem. This is followed by the development of the moment system of
equations through Ito’s formula. We then formulate the Gauss-Galerkin method
and the power series method. Also a summary of other previous methods is given
and an improved numerical algorithm is included. A detailed study of the power
series method is presented and some convergence theorems are established for the
finite interval case. The main convergence result (Chapter 3, Theorem 3.1) gives a
uniform approximation for any j-th moment of the process. The other convergence
result (Chapter 3, Lemma 3.4) gives a characterization of the construction of the

approximating functions.

In Chapter 4, we present the various genetic models for which computations have
been carried out. Several numerical algorithms for the Gauss-Galerkin method and
the power series method are presented. We include a detailed comparison between
the Gauss-Galerkin method and the power series method, for these models. We then
give a steady state analysis and a discussion concerning our numerical results. We
end Chapter 4 by providing concluding remarks and an outline of some suggestions
for future work. An open problem is the approximation of joint distributions of
the process X. A more difficult problem is the approximation of the transition
probability function. An examination of the conditional moments, which is done

at the end of Chapter 4, may lead to some useful ideas for the approximation of

5

the transition density function (when it exists). The Appendix contains numerical
simulations, results related to several examples, and the boundary classification.
According to Feller, the boundary behavior of a diffusion process can be classified
in two categories: accessible and inaccessible. Furthermore, an accessible boundary
can be either an exit or a regular boundary. An inaccessible boundary can be either
a natural or an entrance boundary. Discussions of such classifications as how they

relate to our work are given in the Appendix.

Chapter 2

Statement of the Problem

In this chapter, we introduce some definitions and notations needed for the for-
mulation of the problem. Some known results from both stochastic analysis and
numerical analysis are included. The stochastic analysis results can be found in
Ethier and Kurtz, Wentzell and Arnold. The numerical analysis results can be
found in Stoer and Bulirsh, Stroud and Krylov. Then we present the statement of

the problem.

2.1 Background Material from Probability

Let X = (Xht Z 0) be a stochastic process defined on a probability space (9, f, P)
with state space E(E C R), and let f}=0{X,,s S t}. Let B(E) be the Borel a-
algebra of subsets of E. Then X is a Markov process if the following property holds
(almost surely)

P(X,+,|}}) = P(X,+, 6 I‘lX,),VI‘ E B(E),Vs,t Z 0. (2.1)
Note that 2.1 is equivalent to the following
E[f(Xt+.)I-7::] = Elf(Xt+a)lth,Vf 6 B(E),V8,t Z 0, (2-2)

where B(E) is the space of all bounded measurable functions on E equipped with
the supremum norm. Equation 2.1 is equivalent to saying that ”the past and the

future are statistically independent when the present is known”.

6

7

A function P(.,.,.) defined on [0, 00) x E x B(E) is a time homogeneous transition
function if

1. P(t,x,.) is a probability measure on B(E), Vs,t 2 0,Va: E E,

2. P(O,x,.)=6, (Dirac’s delta function at x), Va: 6 E,

3. P(., .,I‘) e B([0,oo) x E),VI‘ E B(E),

4. P(t + s,x,I‘) = [E P(s,y,I‘)P(t,x,dy),VI‘ e B(E),Vs,t _>_ 0,Vx 6 E.

A transition function P(.,.,.) is a transition function for a time homogeneous

Markov process X if
P(XH-alft) = P(S’Xh F),V3,t 2 09w 6 B(E)-

Let X be a time homogeneous Markov process with transition function P(.,.,.).

Then X will be called a diffusion process if X has continuous sample paths and if O
Vs 2 0,1: 6 E, e > 0, the following properties are true

5. lim,.., 5-37) f|z_y'>€(:r — y)P(t - 3,1, dy) = 0,

6. lim,_., 5:15;)- flz—yISe(x — y)P(t — s,a:,dy) = a(:r, s),

7. lim,_., “—15 flz_ylzse($ — y)P(t - s, x, dy) = b2(:r),

for some measurable functions a and b. The function a is called the ”drift” and
the function h is called the ”diffusion”.

2

The standard Brownian motion with mean 0 and variance 0 is a simple example

of a difl’usion process with drift 0 = O and diffusion b = a.

For a given time homogeneous transition function P(., ., .) and for every f e B(E),

we can define a family of operators (Tht 2 0) as follows

(inmxx) = (mama)

the family (Tht Z 0) defines a measurable contraction semigroup on E (by the
Chapman-Kolmogorov property). To each diffusion process with coefficients a and
b is assigned the second order differential operator

d 1 d2

= — e—2—
L_adx+2b dx2’

8

L( f) can be formally written for every twice differentiable function f and is deter-
mined by a and b. L is known as the Kolmogorov backward operator. Its formal

adjoint operator L', is known as the Kolmogorov forward operator.

The Kolmogorov Backward and Forward equations

 

Let X be a time homogeneous diffusion process with state space E, generated by
the stochastic differential equation dX, = a(X,)dt + b(X,)dB,, X0 being given.
Let a and b be smooth enough and let p(t,:c, y) denote the probability transition
density of X t. Then the right hand side of relations 5 and 6 above coincide with the
coefficients a and b involved in the stochastic differential equation above. Moreover,

the density p(t, x, y) satisfies the backward equation

W = (L).p(t.x.y). (2.3)
p(0.w,y) = 6r.-.» (2.4)

with suitable boundary conditions on p(t, 2:, y). The density p(t, 2:, 3]) also satisfies
the forward equation (Fokker-Planck equation)

9395—5439 = (L‘),p(t,x,y), (2.5)
p(0,$,y) = 6(2—y): (2.6)

with appropriate boundary conditions on p(t, 2:, y).

Equations 2.3, 2.4 plus boundary conditions and 2.5, 2.6 plus boundary conditions
form two parabolic initial boundary value problem. Since the process X may ex-
hibit various types of behavior at the boundaries of E, it is necessary to impose
boundary conditions on the density p(t, x, y). For a more detailed discussion of the

appropriate boundary conditions, see Chapter 3 and the Appendix .

9
2.2 Background Material from Numerical Anal-
ysis

We now present some results from numerical analysis. These results can be found

in Stoer and Bulirsh, Stroud, and Krylov.

Lemma 2.1 Let ,u be a probability measure on ([0,1], 8) such that its support does
not reduce to a finite set. Then, for every n (positive integer) there exist nodes

{3”, 1 _<_ k 5 n} and weights {whml _<_ k g n} with the following properties

1. The nodes are unique and lie in (0,1). They are the zeros of p" (where p“
is the n — th orthogonal polynomial associated with the measure u, described
below).

2. The weights are unique and positive.

3. V f 6 C2”[0, 1], fol f(:c)dp(a:) = 22:1 wk,nf(xlc,n) + €n(f),

where on is an error term given by

f‘””(£)

1 2
e..</)= (2n), /0 one» due).

The orthogonal polynomials associated with the measure u are uniquely determined

by the following three term recursion,

PM = (17 - 5i+llpi - 7511791,

p0 = 1, where
($11.31).) .
6,-+1 = ———,i=0,1,...,
(pivpti
2 (piapi> °
7.- = -———.2=0.1,---,
+1 (Pi—1,1?i—1)
’71 = 0-

The measure p..(.) = 22:1 wk,,,6k,n(.) is called the Gauss-Christofi'el measure asso—

ciated with the measure p.

10

Lemma 2.2 Let the hypotheses of Lemma 2.1 be satisfied. If we define the sym-
metric tridiagonal matrix Jn = Dn + Un + (Un)", where D” and U” are given by

 

 

 

 

(610 o o)
0.520 o
Du: I: I .
00...6,,_10
(00... o 6,.)
Manama)
0073...0

U“: ...
00...0'y,,
(oo...o 0)

Then it is known (for a proof, see Stoer [4]) that Jn has n distinct eigenvalues which
are the zeros of P" (hence the nodes of the measure an). Moreover, if we let A.-
be the i — th eigenvalue of the matrix J,, and v,- be the corresponding normalized

eigenvector, then the weight 212,-," is given by
10.3.. = ((0002.

where (v,)1 is the first coordinate of the eigenvector v,-.

2.3 Statement of the Problem

We are concerned in this dissertation with numerical solutions of the probability

law of a class of stochastic differential equations of the type

dXt = 0(Xt)dt+b(Xt)dB¢, (2.7)
X0 = given, (2.8)

where B. is the standard Brownian motion, and the functions a and b are known

as the ”drift” and the ”diffusion” coefficients, respectively. Equations 2.7 and

11

2.8 model many problems in the physical, engineering and biological sciences (see
Karlin and Taylor, Chapter 15). Under suitable assumptions on the functions a and
b, the existence and uniqueness of the solution of the above initial value problem

is known and the state space can be one of the following:

...:

. (-00, +00),

M

. [0, +00),
. [1, r),

4. (l, r),

09

5. [l, r].

where -00 < l < r < 00. In this thesis, we are concerned with the case where
the coefficients a and b2 are polynomials (in the spatial variable 2:). Under more
stringent conditions on a and b (see Chapter 3, section 3.1), it can be shown that the
state space for the diffusion X can be taken to be the finite interval [1, r]. Without
loss of generality, the state space can be taken to be [0, 1]. By Ito’s formula (see
Chapter 3, section 3.1), it can be shown that the following basic equation holds

dE[¢(X.)l
dt

= E[(L¢)(Xr)l. W5 E 62W. 11- (2-9)
D. Dawson suggested a Galerkin type method for approximating the law of X.
by discrete measures. This method transforms the problem of solving the Fokker-
Planck equation for the density of the law of X 1 (when it exists) to one for a system
of nonlinear ordinary differential equations for the ”nodes” and ”weights” of the
discrete measures (see Chapter 3, system 3.15). This method, when applicable,
works for a more general class of models which includes models where a and b2 are
polynomials. In this latter case, system 3.15 can be viewed as a linear system of 2n
equations for the moments of the atomic measure, where n is the number of nodes
of the discrete measure. This linear system of equations may not be closed in the
sense that these equations may involve the 2n - th or higher moments. Lemma

3.1 of Chapter 3 provides a procedure called the Hankel truncation to close this

system. Then the resulting system (see Chapter 3, system 3.17) is a nonlinear

12

system for the first 2n moments of the discrete measure. Systems 3.15 and 3.17
are both nonlinear systems. System 3.15 is more general than system 3.17 in the
sense when the coefficients a and (or) b2 are not polynomials, it is still meaningful
while system 3.17 is not. It should be pointed out that when system 3.17 can be
formulated, the problem of existence and uniqueness of its solution is equivalent to
the problem of existence and uniqueness of the solution of system 3.15, although
the moment system 3.17 is more tractable (see Chapter 3, section 3.2.3). When
system 3.15 has a solution on some t interval of the form [0,T] and T is a positive
constant independent of n, the Gauss-Galerkin method provides good approxima-
tions (see HajJ afar). Unfortunately, there are models in genetics (see models 1 and
3, in Chapter 4) where T is a function of n causing the solution of system 3.15 to
exist only locally. In the latter case the existence of T as a positive constant must
be established numerically. Moreover, in the Gauss-Galerkin method the initial
distribution is required to be nondegenerate (for more details, see Chapter 3 of the

thesis).

Our main aim in this dissertation is to present a new method, the power series
method, for which theoretical and numerical studies are made (see Chapter 3 of
the thesis). The boundary behavior of the process is especially of importance
for the complete formulation of the system of equations for the moments of the
process which is derived based on equation ?. In Chapter 3 we show how equation
3.4 is more basic than the Fokker-Planck equation, since 3.4 does incorporate the
boundary behavior of the process as well as its behavior in the interior of the state
space, while the Fokker-Planck equation governs the behavior of the density only in
the interior of the state space. When specialized to functions of the type ¢(a:) = 1:",
for k E N, provided they are in the domain of L, equation 3.4 yields an infinite
system of ordinary differential equations for the moments of the process. This
system is open since the first n equations may involve n +1 moments or more. The
number of higher moments involved in the equation for any n — th moment depends
on the degree of the polynomials a and b2 (see Chapter 3 for more details). The

approximation suggested by the power series method is simple in nature. It takes

l3

advantage of the power series expansions of the first few moments of the process

X. This method is based on the following observation

L(x") = lsa(:r:).v‘"l + EUC—g-flb2(x)xk'2,

where L is the Kolmogorov backward operator associated with the process X. If
a is the degree of a and B the degree of b2, then L(x") is a polynomial of degree
at most V(a,fl, k)=maa:imum{a + k — l,fl + k — 2}. When the first 11 moments
are analytic functions in t, their respective power series are truncated at the same
index n, i.e. from the infinite series expansion of any i-th moment of the process,
where O _<_ i S V - 1, we retain the first (n + 1) terms. The V finite sums resulting
from this truncation involve (n +1)z/ coefficients. These (n +1)V coefficients can be
found uniquely, by retaining the first (n + 1)V equations of the original infinite open
system. The full details can be found in Chapter 3. Based on the above truncation,
the resulting system has a unique global solution (see Chapter 3). Then we establish
two main convergence results and give a detailed steady state analysis. Our first
convergence result (Lemma 3.4) provides a way of constructing the approximating
functions needed for the power series method. The second result (Theorem 3.1)
concerns the uniform approximation for any It — th moment of the process X by
a sequence of functions. We also provide a more efficient numerical algorithm for
the Gauss-Galerkin method. This numerical algorithm is an improvement over
the existing one which was suggested by Dawson. The latter scheme uses several
inversions of Vandermonde matrices in order to compute the ” nodes” and ” weights”.
As the number of ”nodes” and ”weights” increases, this numerical scheme becomes
ill-conditioned due to the near singular behavior of the Vandermonde matrices.
Our numerical scheme is based on an eigenvalue problem for the ”nodes” and
”weights” and it reduces the ill-conditioning considerably. The proofs of these
theorems require techniques from both analysis and probability theory. Specifically,
the proofs involve the use of differential inequalities and results from the classical
problem of moments. A number of numerical results related to Models 1, 2 and
3 are presented and, whenever possible, compared with the exact solutions. For
Model 1 the exact solution for the density of the law of X , is explicitly known. For

Model 2, an explicit solution for the moments is available. Moreover, for Model 2

14

an explicit solution for the absolutely continuous part of the law of X t is available.
Several computer simulations are also made for all three models. They can be

found in the Appendix of this thesis.

Chapter 3

Numerical Methods

We start this chapter by stating the assumptions needed for the formulation of our
problem. This is followed by a detailed derivation of a system of ordinary differential
equations for the moments of the process X (system 3.9 below). Then a summary
of some existing methods which have been used for the various truncations of this
infinite system of equations is stated. These methods are: the moment neglect
method, the fixed nodes method and the Gauss-Galerkin method. An improved
numerical algorithm for the Gauss-Galerkin method when the coefficients a and b2
are polynomials is then presented. This is followed by a detailed study of the power
series method indicating two convergence theorems. The main convergence result
(Theorem 3.1) gives a uniform approximation (with respect to the t-variable) for
any k-th moment of the process X. The other convergence result (Lemma 3.4) gives

a characterization of the closure problem for the power series method.

3.1 Assumptions and Derivation of the Moments

System of Equations
Consider the stochastic differential equation

dX, = a(x,)dt+b(X,)dB, (3.1)

15

16

where X0 is given, B, is the standard Brownian motion and the coefficients a and
b satisfy the standard conditions for the existence and uniqueness of the solution.
Under the following assumptions, it is known that the solution of equations (3.1)
and (3.2) is a Markov process with state space [0,1] (see Ethier and Kurtz)

A1. a,b 6 C[0,1],

A2. b(:r) > 0,V:r 6 (0,1),

A3. X0 is independent of the o-fields o{B¢;t 2 0},

A4. a(0) 2 O,a(1) 5 0 and b(0) = b(l) = 0,

A5. D(L) contains all the monomials (where D(L) is the domain of the operator
L).

Remark; It is easy to show the assumption b(O) = b(1) = 0 implies that neither
of the boundaries of the state space is regular. Moreover, if in assumption A4 we

have a(0) = a(1) = 0, then it can be shown that assumption A5 holds.

Now we derive the basic equations which will be the starting point of our work.
Let ¢ 6 C 2[0, 1] fl D(L), then by applying Ito’s lemma to equation (3.1) we have

the following equation
1
d¢(Xt) = {a(X.)¢’(X.)) + 552(X:)¢"(X:)}dt + b(X:)dBo

this implies

t 1 t
ax.) - ax.) = /. {a(x..)¢'<x.) + §b2<x.>¢"<x.)}du + /, b<x..)dB.,
by taking the expectation of both sides of the above equation and then applying

Phbini’s theorem we obtain

E[¢(X.)] -E[¢<x.>1 = [E[(L¢)(X.)ldu- (3.3)

Using the strong continuity of the semigroup (T,;t _>_ 0) corresponding to the
diffusion process X, the function: t ——r E [(L¢)(X t )] can be shown to be continuous.

Therefore, by using (3.3), we obtain

fiEMx.) = E[(L¢)(X.)1, (3.4)
X0 = given. (3.5)

17

This equation always includes the boundary contributions, i.e., contributions from
the singular part and the absolutely continuous parts of the law of X,. The deriva-
tion of a similar equation is possible starting from the Fokker-Planck equation. In
this latter case the boundary contribution may not be accounted for. The following

analysis justifies the importance of the boundary conditions.

Consider the parabolic initial boundary value problem 2.3 and 2.4 for the function
p(t,x,y) (see Chapter 2). We multiply equation 2.4 by ¢ 6 D(L), then the resulting
equation is integrated by parts, twice with respect to the variable x. The following

equation is then obtained

%A1p(t,x)¢(x)dx = ((15, L'p) = (1145.1?) + T,

where T is a boundary resulting from integrations by parts. When T=0 we obtain

the system

%A1P(t,$)¢(fir)dx = /01[(L¢)($)]p(t,x)d:r, (3.6)
P(Oflv) = given. (3.7)

These two equations represent counterparts of 3.4 and 3.5. Unlike equations (3.4)
and (3.5), equations (3.6) and (3.7) generally do not include the boundary behavior
of the process X. As an illustration, we give a few examples where T=0 and D(L)

contains all the monomials (for more details see Ethier and Kurtz, page 366).

Brownian Motion: For the standard Brownian motion with mean 0 and variance
t, the corresponding Fokker-Planck equation for the transition density is given by
0100.3. 31) _ 102mm. 31)

6t 2 63/2

pm. :8. y) = 6(x - y)-

,V$,y E (—00,00),Vt 6 (0,00),

It is easy to show that this initial value problem has a unique solution given by

1

1. _ 2
130.1331) = fiexpkgT'z/L)

18

and the law of X. is absolutely continuous with respect to the Lebesgue measure

on (—00, 00). Its density p(t,x) is obtained as follows

p(t, y) = [190.3, y)p(0.$)dx-

Therefore, in this case no boundary conditions are required for the formulation of

the problem.

Other Models: a E 0 and b(:c) = {RH—:5, 7 > 0. When 7 = $3 this model
is one of the Wright-Fisher models which arises in population genetics. For this
class of models, the state space is [0,1] and the boundary term T depends on the
magnitude of the parameter 7 as follows for 7 > £3 T reduces to —%[p(t, 1)¢(1)]. In
this case it can be shown that T 76 0. Therefore boundary conditions at the end
point 1 must be imposed on the function p(t,x,y). For the case where 7 = % one

can show that boundary conditions at both end points 0 and 1 must be imposed

on p(t,x,y)-

In our studies, we consider equations 3.1 and 3.2 with polynomial coefficients whose

state space is the interval [0,1]. Therefore, we set the following conditions

a 19
a(:v) = z: Akx", and b2(:c) = Z Skxk, (3.8)
k=0 k=l

where 0,,6 E N and A055 75 0.

It can be shown when a(0)=a(1)=0, the domain of L contains all monomials (see
Ethier and Kurtz, Chapter 12, page 367). By specializing equation 3.4 to ¢(x) = xi ,
where j 6 N U {0}, we get

dM- t . ° ‘9
—th—(2 = J z Aij+k_1(t) + 91' z: Sij+k_2(t), (3.9)
k=0 k=l
Mk(0) = given, (3.10)

where Mj(t) = E(th) and 0,- = 312;” . The system of equations 3.9 and 3.10 is an

infinite system of linear equations for the hierarchy of the moments of the process

19

X. This system is closed if and only if

osa+j-13j and 0_<_fi+j-2_<_j,VjeNU{0},

this is possible only when a 5 l and ,6 5 2, in which case the exact solution of this

system is given recursively by
t
Mj(t) = [Mj(0) - Mj_1(0) + 91'!) Mj_1(8) exp(—0js)ds] exp(0jt).

If a Z 2 or fl 2 3 then 3.9 is not a closed system and therefore cannot be solved ex-
plicitly. Our task is to close the system 3.9 in such a way that the resulting solution
yields a good approximation to the original system 3.9. Various truncation schemes
have been proposed to close the system. For some models which arise in Genetics,
the derivation of system 3.9 and 3.10 is not possible from the Fokker-Planck equa—
tion. This is an important observation because when the boundary behavior is not
incorporated, the steady state analysis for the process X is not possible. We em-

phasize the importance. of the boundary behavior through the following discussion.

Consider the following class of models which arise in population genetics quite
often:

a(:c) = x’(1 — x)’ and b(:r) = x”(1 — 1:)", where r, 3,2}? and 2q 6 N.

Note that a and b2 can then be written in the following form
a l3

a(:c) = z Akx", and b2(:r) = X 51.3,”,
k=0 k=l

witha=r+sandfl=2(p+q).

For a E 0, we treat the following cases here. The discussion of the more general

case is contained in Appendix.

Case 1: p,q 2 1 (this implies that both 0 and 1 are natural boundaries), we have
T E 0 and
<¢, L'p) = (Lo, p),v¢ e C2[0,1].

20

Case 2: p = %,q 2 1 (this implies that 0 is an exit boundary and 1 is a natural
boundary), T is given by
1
T = —§¢(0)p(ta 0):

which implies

<¢.L*p> = (lap) — §¢(0)p(t.0).\r¢ 6 070,1].

Case 3: p Z l, q = z],- (this implies that 0 is a natural boundary and 1 is an exit
boundary), T is given by
l
T = -§¢(1)p(ta 1):

which implies

(¢.L"p) = sap) — grow, an e C2,, 11.
Case 4: p = q = i; (this implies that both 0 and 1 is are exit boundaries), T is
given by
T = -§[¢(1)p(t. 1) + ¢(0)p(t. 0),

which implies

(45,131?) = <L¢.p> — %[¢(1)p(t.1)+ ¢<o>p(t,o>1,v¢ e 62m. 11.

Cases 2, 3 and 4 require special analysis because of the boundary behavior of the
process. At an exit boundary the density p(t,y) > 0 (see the Appendix) which
implies that T 96 0, for some <15 6 CZ[0,1].

3.2 Summary of Previous Methods

3.2.1 Moment Neglect Method

This is the simplest and most frequently used method. If the state space of the

stochastic process X is the interval [0,1], then all the moments exist (they are all

21

bounded by 1, uniformly in t). The following observation is the basis for this
method: for each fixed t we have

1 1
0 5 M..(t) = E(Xl‘) =/0 x”p(t,:r)dx 5 D(t)]0 as” = D(t);—.1;T —r 0, as n —+ 00,

where D(t) is assumed to be finite and is defined by

D(t) = sup{p(t,:r);a: 6 [0,1]}.

Therefore in system 3.9 it is reasonable to set M..(t) equal to zero for large values
of n. The method consists in setting the (n+1)-th and higher moments equal to
zero. In this way the infinite system of moment equations for the process X reduces
to a closed finite linear system of n equations. There are two basic difficulties with
this approach (or variants of it). The method leads to an unbounded solution and
it does not lead to convergence. For a more detailed discussion, see Chapter 3. As
an illustration of the moment neglect method, we present the following model (this
model will be referred to as model 3, see Chapter 4 for this and two additional

models)

050 and b(:1:)=:r\/1-:r,

for this model system 3.9 yields the following equations

M120) = BkleU) - Mk+l(t)l1Vk E N U {0}.

for k large enough, MI,“ is neglected leading to the following finite system of linear

equations

m;(t) = 0,-[mJ-(t) — mj+1(t)],Vj = 0, 1,2, ..., k — 1,
"15.0) = 9kmk(t).
the last equation implies
mk(t) = mk(0)exp(0kt).

Consequently,

mic—1(t) = gk—llmk—l(t) - mk(t)l. implies

22

t
mic—1U) = [mic—l - aka/0 mk(3) eXP(-0k-13)d3l eXP(9k-1t).
and in general we have the following recursion valid for j = 2, 3, ..., k + 1

mj-1(t) = [mj_.1 — 01--1A:mj(s)exp(—0,-_1s)ds]exp(0j_1t).

This method may give satisfactory results only for very small values of t. It is
observed (see Chapter 3) that the first few moments have exponential growth and

exceed 1 even for small values of t hence giving misleading results.

3.2.2 Fixed Nodes Method

In this method 11 nodes are fixed according to a given distribution (the distribution
of X0). Then 3.9 can be made closed by expressing higher order moments in terms
of lower order moments as follows. If we let 2:1(n,0),:r2(n,0), ...,xn(n,0) be the n
nodes of the Gauss-Christofi'el measure corresponding to X0 (see Chapter 2, Lemma
2.1) and let

”(nat) = Z wi(ntt)635(fl,0)t

i=1
then in system 3.9 when the moments of the process X are replaced by the moments
of the measure u(n,t), it can be shown that the weights of u(n,t) satisfy the

following initial value problem
VW’ = BW, (3.11)
W(n,0) = given, (3.12)
where (W(n,t))T = (w1(n,t),w2(n,t),...,wn(n,t)), and u(n,0) being the Gauss-
Chfistofiel measure associated with X0,
V=v.-, ={z{,.‘=1,2,...,n;j=0,1,...,n— 1},

and
B = {bgj} = {[L(1")](x,-);i = 1,2, ...,n;j = 0,1,...,n -1},

23

this system has a unique solution given by
W(n, t) = W(n, 0) exp[t(V"l)B].

To express the higher moments of p(n, t) in terms of the lower moments, we proceed

as follows. Define the matrix V(k) by

 

 

( (21):" (2:2)“ (33)“ my» \
(2:1)Im+l (x2)kn+l ($3)lm+l . . . (1.")!an
V(k) =
($1)lcn+n—2 (x2)kn+n-2 . . . (xn_l)kn+n—2 (xn)lcn+n—2
\ (1'1)""+."‘l ($2),m+n—l . . . (xn_1)k”+“‘1 (xnyflH-n-l }

Note that V(0) is the usual Vandermonde matrix. We define the following vector

m(k, t) = (mien: mkn-i-l: "'a mkn+n-l)a

where m,-(n,t) = 233:, wj(n,t)[xj(n,t)]. Then for each k=0, 1, 2, ..., we obtain the
following system
m(k,t) = V(k)W(n,t).

For k=0, we have
W(n,t) = V(0)“M(0,t).

which implies
m(k,t) = V(Ic)[V(0)‘1M(0,t)].

The fixed nodes method is easy to implement and often yields satisfactory approx—
imation of the moments of the process X for small values of t. This method is
similar in nature to the well known Newton Cotes method for approximating inte-
grals. Unfortunately, we do not have a convergence theorem of the approximating
moments to the true moments, as n —> 00. In the following section we present the

Gauss-Galerkin method. In this method the following assumption is required:

A6. Equations 3.4 and 3.5 above, uniquely determine a measure for all t E [0,T].

24

3.2.3 Gauss-Galerkin Method

We start by reviewing the Gauss-Galerkin method. Then we shall present a new
numerical algorithm for the implementation of the method. Let um, be the Gauss-
Christoffel measure associated with the law of X, (see Chapter 1). Then system

3.9 can be written as follows
2,1,1; n..<t)(e...(t)>‘1 = [21 ...,(r)L(g,,(.)y] + 6...”,

for all i=0, 1, 2, ..., n, . This system cannot be used to find the weights 17,, J(t)
and the nodes (“J-(t) because of the error term involved in the above system. The
Gauss—Galerkin scheme is obtained by dropping the error term in the equation

above. This way, as stated in Chapter 2, we obtain a sequence of discrete measures

n
“mt = Z wnJ(t)61’u.j(‘)’

1:1
where the weights wad-(t) and the nodes 1:" J(t) are solutions of the nonlinear finite

system of dynamic equations

i123, WWW”)? = gaseous-(0):], (3.13)
wnJ(0) = 77n.j(0) and mud-(0) = 6mm). (3.14)

for all i=0, 1, ..., 2n-1. The measure
11...: = z wn,j(t)6:r.,,-(t)
i=1
is called the Gauss-Galerkin approximation of the law of X,. If the nodes and the
weights are differentiable functions of t, then the dynamic nonlinear finite system

3.13 and 3.14 can be written as follows:
A(n,t)X(n,t)' = D(n,t)X(n,t), (3.15)
X(n, 0) = given, (3.16)

where the initial conditions and matrices A(n,t) and D(n,t) are defined as follows:
X(n1 0) = (nn,l(0)a "713(0): "-1 nu,n(0)a £n,1(0)1 £n,2(0)1 "'3 £n,n(0))a

25

A A
A(n,t) = 11 12 ,
A21 A22
D 0
D(n,t) = 1 ,
D2 0

A11 = {(xnj)i}n>(n:

A12 = {wnJKanl‘l'hxm
A21 = {(xnj)i+n}nxvn

A22 = {wnJ[($n.j)‘+"l'}nxm
DI = {L($‘)($n.i)}nxm
D2 = {L(x‘+")(xn,j)’}..x...

The existence and uniqueness of the solution of system 3.15 has not been settled
except for some special cases. In case the coefficients a and b2 are both polynomials,
system 3.15 describing the dynamic behavior of the nodes and weights can be viewed
as a finite linear system of equations for the moments of the discrete measure. The
latter system is a truncation of the infinite linear system of ordinary differential
equations 3.9 which is satisfied by the true moments of the process X. However, the
system 3.15, when viewed as a system for the moments of the discrete measure, is
not closed in the sense that the equation for the (2n-1)-th moment involves 2n-th
moment or possibly higher moments. A procedure called the Hankel truncation is
used to close this system. The resulting system is a nonlinear system of 2n equations
for the moments of the atomic measure. To do this, the Hankel determinants of
order higher than n are set equal to zero. The Hankel determinants are defined as

follows, for each r 6 N U 0,

mo m1 . . . m,

m1 m2 . . . m,.+1

 

 

mr mr+1 e e 0 m2,-

26

and
m1 771.2 . . . mf+1
m2 m3 . . . mr+2
r, =
mr+1 mr+2 - ° - m2r+l

 

 

Lemma 3.1 Given a sequence of numbers {m,-,i Z 0} there exist a unique atomic

measure whose support consists of exactly 11 nodes if and only if
I. A, > 0,Vr E 0,1,...,n — 1,
2. I‘, > 0,‘v’r E 0,1,...,n — 1, and
3. A, = Pr = 0,Vr 2 11.

Both Dawson and HajJafar assumed (i) and (ii). The assumptions are then ver-
ified numerically (a theoretical proof is not available). As a consequence of these
assumptions the higher moments are expressed in terms of the lower moments as
follows

An(t) = An-1(t)m2,,(t) — g[mo(t), m1(t), ..., m2"-1(t)], where g is a polynomial. By
setting An(t) = 0, we obtain

g[mo(t), m1(t), ..., 7712"_1(t)]
A.._1(t)

 

m2n(t) = = f[m0(t), m1(t), ..., m2"-1(t)].

Let’s give an illustration of the method through model 3 (i.e., a E 0, b(x) .-=
le - x. Then the system of equations for the moments of the process takes the

following form

m;(t) = 0,-[mj(t) - mj+1(t)],Vj = 0,1, ...,2n - 2,

man—1(t) = 62n-llm2n-l(t) - flm0(t)r ml(t)r '"r m2n-l(t)]a
which can be written as

m'(n,t) = A(n)m(n,t)+F(n,t), (3.17)
m(n,0) = given, (3.18)

27

where
m(n,t) = (mo(t),m1(t), ...,m2n_1(t)),
F(n,t) = (0,0, ...,0,f[mo(t),m1(t), ...,m2n_1(t)]),
.401) = {055},i,j = 0,1, ...,2n — I,
where
0,- if i=j
aij = -0,~ ifj=i+1

0 otherwise.

Lemma 3.2 The system 3.17 and 3.18 has a unique solution if and only if there
exist a neighborhood [0, T] such that Vr = 0,1, ..., n - 1 and t E [0, T], the following
holds

1. A.(t) > 0,
2. I‘,(t) > 0.
(The proof is straightforward).

Lemma 3.3 Let the 2n functions {m,-(t),0 5 j 5 2n — 1} obtained in Lemma
3.1 be the first 2n moments of an atomic probability measure p(n,t). Then Vt E
[0,T],V1,ZI 6 C[0,1], then as n —r 00, the following is true

[01 mm t) —> [01¢(X.)dP.

The nodes and weights of the discrete measure p(n,t) can be found through the
usual inversion of the well known Vandermonde matrix. This matrix is given by
V(0) of section 3.2.2 above. It is easy to see that the determinant of V(0) is directly
proportional to

n

H [(121.40 - x..,,-(t)]2.

5J=1.i¢j
This quantity may tend to zero as It tends to infinity. As a consequence V(0) may

become nearly singular, causing the problem of solving for the nodes and weights

of the measure u(n,t) to become ill-conditioned. An alternative method which

28

considerably reduces the ill-conditioning when computing the nodes and weights is

the content of the numerical algorithm below.

An Improved Numerical Algorithm for the Gauss-Galerkin Method

 

Our numerical algorithm for finding the nodes and weights of the Gauss-Galerkin
measure p(n,t) is based on Lemma 2.1 and Lemma 2.2 of Chapter 2. Lemma 2.1
gives a relationship between the orthogonal polynomials associated with p(n, t) and
its first 2n moments {mJ-(t), 0 5 j 5 2n - 1}. Lemma 2.2 provides a way of com-
puting the entries of the tridiagonal matrix Jn. We shall now present the detailed

derivation of the algorithm.

Given the 2n moments {m, (t), 0 5 j 5 2n — 1}, the inner products involved in the
construction of the matrix J" can be expressed in terms of the given moments as

follows

(P01P0> = m0)
($130,170) = ml,
in general, given

2.‘ ,
(Path) = Zagz')m2,-_j,where0 5 i 5 n — 1,
i=0

and the a?” are known. Then

2.‘
(xpirpi) = Z a§2i)m2i-j+l-
j=0
Through the three term recursion formula we also have the following
($191.19.) - 5i+1<Pi.Pi) = 0 =>
(Pi+1.P.'+1) = aligning + - -- + Gigantic - 25i+1(082i)m2i+1 + - - - + Ohiomll +
(6.,1)2(a32‘)m2. + . . . + agimo) —

(Palms-1) + (’Y.‘+1)4(Pi—1.Pi-1)-

29

We now describe a one step computational algorithm from t = 0 to At. This
computational algorithm is used to compute the approximate values of the first 2n
moments of the process (the process arising from model 3) at time At using the
initial conditions {Mk(0),0 5 k 5 2n - l}.

1. Given {Mk(0),0 5 k 5 2n — 1} we have (for the specific example under
consideration)
M2“) = 0..[M,.(t) - Mk+1(t)l

which implies

MUM) = Mk(0) + othtm) - Mk+1(0)l-

2. With the new mk(At), k = 0,1, ...,2n — 1, we compute 6,,(At), ..., 6,,(At) and
7,,(At), ..., 7,,(At) through Lemma 2.2 of Chapter 2. Then we form the matrix
J..(At). The eigenvalues of Jn(At) and their corresponding eigenvectors are
then found. This yields the nodes {x,-(n, At), j = 1, 2, . . . , n} and the weights
{w,-(n,At),j =1,2,...,n}.

3. Using the first iterates, we obtain the second and higher iterates from {mk(At), k =

0,1,. . . , 2n — 1} through the following formula
mk((i + I)At) = mk(iAt) + 0k(At)[mk(iAt) - mk+1(iAt)],
where i 2 1.

This numerical algorithm has proved to perform better than the algorithm proposed
by HajJ afar and Dawson in which several matrix inversions are required causing the
problem of finding the weights and nodes to become ill-conditioned. Our algorithm
does not require any matrix inversion. Therefore, our numerical algorithm is more

suitable when higher values of n need to be considered.

3.3 The Power Series Method

When a and b2 are polynomial functions, system 3.9 of equations for the moments

of the process X is closed if and only if \

30

OSa+j-1_<_j and 05fl+j-2SJ.VJ'€NU{0}.

this is possible only when a 5 1 and fl 5 2. In the latter case the unique solution
of 3.9 is easily obtained. From 3.9 it is clear that Vi E N, M.- 6 C°°[0,1],VT > 0. In
the Gauss-Galerkin approach the existence and uniqueness of the truncation of 3.9
is known to hold only on a t interval of the form [0,T(n)], where T(n) is a positive
function of n. For more details we refer to Chapter 4. The power series method is
based on a different kind of approximation. The two approximations (the Gauss-
Galerkin and the power series) are not in the same sense. The Gauss-Galerkin
method, when applicable, yields an approximation in the sense of convergence in
measure. The power series method gives an approximation in the sense of uniform
convergence and is applicable only when the coefficients a and b2 are polynomials.
For a more detailed comparison of the methods, we refer to Chapter 4 where several

examples are illustrated.

If a and fl are integers such that a 2 2 or fl 2 3 then 3.9 is not a closed system and
therefore cannot be solved explicitly. Our task is to close the system 3.9 in such a
way that the resulting solution yields a good approximation. For each i 6 N U 0,
the differential equation involved in system 3.9 is

4M“)

a [3
dt = j Z Aij+k_1(t) + 0]: 2 Sij+k_2(t). (3.19)
le=0

Ir=1

When a 2 2 or fl 2 3, system 3.19 is not closed. It is impossible to explicitly solve
equations 3.19. It is assumed that all the moments have power series expansions
in t. In fact, as we shall show below, it suffices to assume the first few moments to
be analytic in t. We then approximate these moments by polynomials. The details

are given below.

Depending on the parameters a and ,6, the following cases are considered:
(i) ,6 at a +1, (ii) ,8 = a +1. In case (i) there are two possibilities, either fl > a +1
or fl < a + 1. For simplicity, we shall treat the situation where fl > a + 1. The

case where 3 < 01 + 1 can be treated similarly. We now consider case (i).

31

Equations 3.19 above can be written as follows:

. 544-2
flit“) = _Z Bs.M(t) (3.20)

 

where Bid = B,-J(i,0.-,A,-, 5,) is given by
B _ iAj...'+1+ 0,512.5.” If! -' 1 _<_j S a +1— I
h] 0;Sj-5+2 ifa+i5j5fl+i—2.

Suppose that the moments M1, M2, ..., M54 are analytic functions on [0,T]. Then

for i=1 equation 3.20 gives

All“ )' 21—2 BIJMJ'
Bl’p- ’

 

Mrs-Kt):

which implies M34 is an analytic function. Note that 81,5..1 ¢ 0. In fact B{,fi+i_2 =
0,3,3 9f 0, Vi E N. For i=2, equation 3.20 leads to

M50) - 2f;1132,ij

MW) = 32x .

 

which implies Mp is an analytic function (since M54 is analytic). More generally,
equation 3.20 implies M,- is an analytic function for all i 2 B — 1. Therefore the

following expansions are possible

M,(t) = Ecj(i)tj,Vi = 1,2,3, ...,n - 2,
j=0

where j!c,-(i) = Mij)(0), for allj 2 0.

Polynomial Approximation Fix any integer n and define the following polyno-

 

mials
n . tj
m,(n,t) = 29(1):, (3.21)
i=0 3'

where i = 1, 2, ..., fl — 2. We then consider the initial value problem

d , t 5+"?
"'(n )= 2.: B,,,m,-(t) (3.22)
j=i- -l

32

subject to the condition
"11(0) = M'(0)1

for all i = 1, 2, ..., u, where u = n(fl - 2). The number of unknowns cj(i) involved in
power is (n +1)(fl - 2) (note that we have co(i) = M,(0), for all i = 1,2, . . . ,fi - 2.
System 3.22 is not a closed system. In the following lemma we suggest a way to

close it.

Lemma 3.4 Using the system 3.22 the coeflicients in (3.21) can be determined by
the the initial conditions M,-(0), i = 1,2, . . . , (n + 1)(fl — 2).

3:991; The coefficients involved in 3.21 are

61(1) 61(2) 61(3—2)
(32(1) 62(2) 02(3—2)

cn(1) c,,(2) cum—2).

The coefficients co(1), co(2), . . . , Com — 2) are uniquely determined since by defini-
tion we have co(i) = M.-(0), for all i = 1, 2,... ,fl — 2. To determine the coefficients
c1(1),c1(2), . . . ,c1(fl — 2) we consider the first fl — 2 equations involved in sys-
tem 3.22. These equations are then evaluated at t = 0 to yield the following

fl-i-i-2

61(1) = Z B.ij(n,0),

j=i-l

for all i = 1,2, . . . ,fl - 2. To determine the coefficients c2(1),c2(2), . . . ,c2(,6 - 2)
the first (fl — 2) as well as the next (,6 — 2) equations are needed. This is clear from
the following equation which is obtained by differentiating the i — th equation in

system 3.22

c12m__,-__(n,t)= “"2 dm (n, t)
cit—T 2 B‘”— —_J__d ’

j=i- l
which can be written as

13+i-2 [NJ-2

2 Bid 2 Bjkmk(n,),t

j=i-—l k=j- -1

33

then “1,4

2lc2(i) = X j = i — 13+‘-2B.- J 2 B-,,.M,.(0).
k=j-l
Continuing in this fashion, all the coefficients {c,-(i),1 5 j 5 n, 1 5 i 5 fl — 2} are

uniquely determined through the relation

 

m2+2_5(n,t) - l:=li+2-B Bi+2—fl,lcmk(na t)
Bin—M
The total number of equations required is n(fl - 2). The number of initial condi-

tions used is (n + 1)(fl - 2). This completes the proof.

mi(nrt) =

Now we turn to the question of approximation of the moments of the process X.
The following Theorem provides a uniform approximation of any I — th moment of

the process X.

Theorem 3.1 Suppose assumptions A1 through A5 hold and assume that the first
fl - l moments of the process X are analytic. Fix any integer l and consider the
sequence of functions {m,~(n,t), n _>_ 55—2}. Then for any To 6 (0,T), this sequence

converges uniformly (in t) to M;(t), on [0,To], as n -r 00.

Emf; The following elementary result from the theory of power series (see Smirnov,
page 150) is needed. This result, when applied to the analytic moments {M,~(t), 1 5
i5 fl—2}, states the following for each i = 1, 2, . . . ,fl—2, VI 6 NUO, VTo E (0,T),

HM! - mi(n.t)ll -: 0." -* 00 (I)

where [I I] is the supremum norm on [0,T].

Let I 6 N, then we have

1. If 1 5 l 5 fl — 2, then by the definition of the functions m,(n, t), the theorem
holds by applying relation (I).

2. If I _>_ fl — 1, then the (l — (B — 2)) - th equation from system 3.22 above gives

1 l—l

mi—a+2(n,t)= Z BI—a+2.jmj(",t)= Z Bl—B+2.jmj(nvtl+Bl—B+2Jml(”at)
j:i-B+2 j:l—fi+2

34
which implies

mi— s+2(n t) 2j=l—fl+2 B1_,9+2ij(n t)
Bl—fi+2.l

 

m1(n,t)- —

Consider the first l— 3 + 2 equations from system 3.22. The first of these equations

leads to :fi_2
m’1(n,t)— oBlJmJ-(n, t)

31).-

 

mfi-1(n1t) =
which implies

m’ n ’ - 2 1J mJ-n, j
lMa—1(t)-ma—1(t)l=|[ 1”) Mm] 2;? l ( t) Mun,

 

therefore

 

 

m’ n,t M’( t()l m n,t M-t
1( B)- l l‘l'lZlBldl 1(B)- J()[—r0,n-+OO,
lfi-l Bl,B-l

9

[M5_1(t)—m3_1(t)| S I

by part 1 above.

From the second equation (of the first I — fl + 2 equations) we have

m2(nrt) - Egg-:11 B2ij(nrt) =
32,8

 

ma(n.t) =

m’2(n,t) — 233:21B2ij(n, t) _m’1(n, t)— 2%: BlJm,(n, t)

 

 

Again because the moments of the process satisfy system (3.9) we can write

m2(nrt) _ M2(t)|+
B2 ,.

 

WW) - W(tll S l

5‘2 m-(n t) - M-(t)
B . J i J +
12:; l 2d B2,fl l

 

,m’1(n.t) - M10)
B2’fllBl’fi_l

$31M, 0- M.-<t)l _,0 .._. 00
j=0 ‘1 B2rflBlrfi-l

 

 

35

by part 1.

Using the same procedure we can write in general
HM; — m,“ —i 0,12 —> 00.

This completes the proof. The case where [i = a +1 requires the following assump-
tion

iAa + 0‘55 ¢ 0,Vi E N U {0},

and a convergence theorem like Theorem 3.1 can be proved in a similar way.

Now we give sufficient conditions for analyticity of the moments of the process X

and provide some examples where the moments are analytic.

Lemma 3.5 Suppose the law of X, is absolutely continuous with respect to the
Lebesgue measure on [0,1] for allt E [0,T] and let {((bi, A;),i 6 N} be the eigenpairs
corresponding to the differential operator L‘. Moreover, assume p(t, y) the density
of X, has the expansion Egon-(Md where i!p.~(y) = Who, and the initial
density p(0,y) lies in the span of finitely many eigenfunctions of L‘. Then the

following are true
1. The formal series above for p(t, y) converges everywhere.
2. The [C - th moment of the process X is an analytic function on [0,T].

£1991; Consider the Fokker-Planck equation for the density p(t, y),

ap(t’ y) _ t
at '— L p(t, y),

where p(O, y) is given. Then from the expansion for p(t, y) and the Fokker-Planck

equation we must have the following relation

i ... (i)
Pi(y) = (zyllmo = (L) if(0,y),i = 0,1,...

 

Since by our assumption p(O, y) = 2:520 AJ-qSJ-(y), then for all i = 0,1, 2, . . ., we get

= ((L')<z°5)>:§;oAj¢j<y)

i!

 

Pd?!)

36

which implies

 

21:0 A,(/\,-)‘¢,-(y)
1).-(y) = J ,,
(since L‘dy = qubj).
Therefore
ling/ltd _<- 2: lpi(y)t'l— — 20% laAJ(AJ) ¢j(y)l < NmfxllAjll¢jllwl :(lANlt)
i- — j=0 i=0
(where IAN] 2 |/\j|, for j = 0,1,.. ., N. Since the eigenfunctions qSJ- are continuous

on [0,1], we have
N N N
guanine} s rgggc{IA.-I}{r;nggcn¢jno.} = c
where C is a finite constant. Therefore

2 lpi(y)lt' < 02— (——._|AN!lt ) — "C €Xp(|ANIt) < 00, Vt > 0.

i=0 i=0
Hence the density p(t, y) is an analytic function of the variable t. Now we shall

prove part 2) of the Lemma. Let k E N, then

Mkt()=/p(ty)y"dy=/219431)? "=Z/oy"dy,

i=0

by the Legesgue dominated convergence Theorem. Moreover,
I [01p1(y)y"dyli _<_ {/01 lpe(y)ldy}i,
which implies
JlggfiupU/y My )yikdz/I} < 11m sup{/ Ips(y)|dy} < {/0 (RWY-ids} =R+€

for all 1 large, where R 18 the radius of convergence for the series corresponding to
the density p(t, y) and e is any arbitrary positive number. Thus by letting e —+ 00,

we can conclude that Mk(t) is analytic. This completes the proof.

Examples of Models which possess Analytic Moments A few examples where

 

all the moments of the process X involved in the model are analytic functions are

as follows:

37

Example 1:
a(a:) E 0 and b(x) = ‘/z(1 — :c),Vx 6 [0,1].

It is easy to see that system 3.9 for the moments yields

Mi“) = 0,-[M.-_1(t)—M,-(t)],
M;(O) = given.

This system can easily be solved explicitly. The solution has the form
t
Mm) = {[M,(0) — M._1(0)] + 9,. [0 141-1(3) exp(a,,s)} exp(—6kt).

The solution clearly shows that all the moments are analytic in t.

Example 2:
It is equally easy to show that the moments of the process involved in model 4.1 are
all analytic functions in t. An explicit solution for the moments system of equations

is available.

Chapter 4
Some Models From Genetics

In this Chapter we present the various genetic models for which numerical compu-
tations have been carried out. Several numerical algorithms for the Gauss-Galerkin
method and the power series method are presented. For each model considered,
we solve the truncated system of moment equations (system 3.9, Chapter 3), as
suggested by the power series method. This is followed by a detailed analysis of
the steady state of the process through the power series method. We then present
a detailed comparison between the Gauss-Galerkin method and the power series
method. We end Chapter 4 with a discussion concerning the numerical results
and future work. These numerical results can be found in the Appendix. In the
following section we present a brief description of the various genetic models under
consideration. For a more detailed description of how these models are derived we

refer to Ewens, Karlin and Taylor, and Kimura and Crow.

4.1 Description of the Various Models

The following models where b(x) = Ax”(1 - x)? are considered:

1. ModellzaEOandp=q=A=L

2. Model2: aEOandp=%,q=%,/\ l

3. Model3: a50andp=1,q=%,/\ 1

38

39

4. Model 4: a(:1:)= sa:(1 — :r)[h + (1 — 2h):c] — 112: + v(1— 55),]; = q = (,\)2 =%

where s, h, u, and v are parameters.

We shall consider two special cases of Model 4. Model 4.1 corresponds to the case
where s = 0, h = 0.5, u = 0.5, and v = 0.5. Model 4.2 corresponds to the case
where s = 4.8, h = 1.3, u = 0, and v = 0. Our numerical results concerning these
two models will be compared to the numerical results obtained by J. Keller and R.
Voronka.

The four models above are important stochastic processes which describe the change
of gene frequency in a Mendelian population. These models are of interest in pop-
ulation genetics. The models describe the behavior of the gene frequency of a
given population, as it fluctuates in time. Consider a pair of alleles A and a with
respective frequencies 3 and 1 — 2:. Then Model 1 describes the case of random
fluctuation of selection intensities. Model 2 describes the simplest situation. It
describes the case which is known in population genetics as the process of random
genetic drift due to random sampling of gametes. Model 3 is a model of mathe-
matical importance. In this model, the moments system of equations is not closed.
Model 4 includes the effects of selection, mutation and migration as well as random
drift due to the finite size of the population (note that through a time change, the
population size can be considered to be equal to 1). In this model the parameter 3
denotes the selection advantage of the zygote AA over the zygote Aa. The param-
eter h is such that the quantity sh denotes the advantage of the zygote Aa over
the zygote aa. Moreover, the parameter u represents the effective mutation rate of
A to a and the parameter v represents the effective mutation rate of a to A (the
effective mutation rates include the effects of both mutation and migration). In
our considerations, Model 4.1 (which corresponds to s = 0, h = 0.5, v = 0.5, and
v = 0.5), and Model 4.2 (which corresponds to s = 4.8, h = 1.3, u = 0, and v = 0)
are well known models in population genetics. For more details see Ewens, Kimura
and Crow, and J. Keller and R. Voronka. We note that all these models satisfy our
assumptions. In the next sections of this Chapter, we shall describe the numerical

results of the theory presented in Chapter 3 as it relates to the four models above.

40

4.2 Solution of the System of Moments

In this section we illustrate the various numerical methods we considered earlier.
We shall apply these numerical methods to Models 1, 2, 3, and 4.1 and 4.2. This
is followed by a discussion of the relative merits of these methods as suggested by
the numerical results. The numerical methods considered are the moment neglect

method, the Gauss-Galerkin method, and the power series method.

4.2.1 Moment Neglect Method

The moments system of equations for the four models are as follows

Modem;

M1,:(t) = ohm/11:0)’2Mk+1(t)+Mk+2(t)l, (4-1)
Mk(0) = given. (4.2)

M93121};
ML“) = 9kle-1(t) — Mk(t)]. (4-3)
Mk(0) = given. (4.4)

Model 3:
Ml’c(t) = 9kle(t)‘Mk+1(t)l, (4-5)
M1,,(0) = given. (4.6)

MQQM

M},(t) = k{(% + v)M,,_1(t) + (£0: — 1) — v — u + sh)Mk(t) + (4.7)
8(1 - 3h)Mk+1(t) - 3(1 - 2h)Mk+2(-t), (4.8)
Mk(0) = given. (4.9)

As an illustration, we apply the moment neglect method (for more details, see

Chapter 3) to Model 3. The resulting system of equations has the following form

m5“) = 92lm2(t)-ma(t)l,

41
m3“) = 93[m3(t) - 77140)],
mix—1U) = on-llmn-l(t) - mn(t)la
mn(t) E 0,
this system can be written in matrix form as follows

m'(n, t) = Nnm(n, t),

 

 

m(n,0) = given,
where
(02—02 0 0... 0)
0 0 —0 0 0
N": o o3 3 o o o
\0 o woos.-.)
and
( m2(t) )
t
m(n,t): "'3,” ,

 

 

l "in—10) /

the unique solution is given by
m(n,t) = m(n,0) exp(tN,,).

When the initial distribution of the process X is uniform on (0, 1), the initial con-

dition is given by

#IH wlh‘

m(n, 0) =

 

 

If
3'... 00.
\

For n=20, the approximations for M2,M3,...,M8, are summarized in the Ap-

pendix. Since the true moments are bounded by 1, it can be seen from the graphs

42

in the Appendix that the moment neglect method leads to misleading results even
for small values of t. When applied to models 1, 2 and 4 the moment neglect method
also gives misleading results. Higher values of 11 gave similar results. Based on this
observation, we do not recommend the use of the moment neglect method or any

of its variants.

4.2.2 Gauss-Galerkin Method

With the initial density being uniform, and with n=5 we compute the initial
nodes $5,1(0), 353(0), . . . , 355(0) and weights w5,1(0), w5,2(0), . . . , w5,5(0) through
the tridiagonal matrix method described in Chapter 3. With these initial values
we solve the moment problem using our improved algorithm which is described in
Chapter 3. The results are summarized in Tables g.2, g.3, g.4.1, and g.4.2, and
Figures 1.1, 1.2, 1.3, and 1.4. Our computation results show that the power series
method and the Gauss-Galerkin method are in close agreement. For models 4.1
and 4.2, both methods support the results obtained by J. Keller and R. Voronka.
However, the method used by J. Keller and R. Voronka is valid only for small values
of t. The Gauss-Galerkin method and the power series method have the advantage
of making the steady state analysis possible.

4.2.3 Power Series Method

In this section we describe the numerical results of the power series method. The
actual numerical results can be found in the Appendix. We have proved the fol-
lowing result (see Theorem 3.1 of Chapter 3), VT > 0, and We 6 N,

lim sup IMk(t) — mk(n,t)| = 0.
"“°° 05:57

Therefore, given any k E N we can approximate (uniformly) the moment M], by
the polynomial mk(n, .) on [0,T], for large values of n. We note that when n gets
larger, more and more initial conditions are needed. Since mk(n, .) is a polynomial
in t, no matter how large n is , we have mk(n,t) -» 00 as t -> 00. Therefore

to analyze the steady state solution, a single polynomial of degree n (no matter

43

how large n is) will not provide useful information when t -> 00. To illustrate
the approximation resulting from the power series method for finite values of t,
we choose the initial distribution to be uniform and fix n = 7. Using the power
series method described earlier, the following results are derived for model 3 (see

Appendix for results concerning all three models),

7
m2(7,t) = 20¢...
i=0

the coefficients c,~, for i = 0,1,..., 7 are found to be

— — 1 — .1. — l.
CO - 313361 "' '13-’62 "' _30903 - 729
_ _a - _2§_ _ __1_9__ _ .2231;
c4 - 315,65 " 7200’C6 — 10800’67 — 2116800'
The rest of the functions m3(7,t), m4(7,t), . . ., m9(7,t) are easily obtained from

the function m2(7,t). They are tabulated in the Appendix. This results are very

close to the results provided by the Gauss-Galerkin method as shown below.

44

4.2.4 Comparison of the Power Series Method and the
Gauss-Galerkin Method

In this section we give a detailed comparison between the power series method and
the Gauss-Galerkin method. We present advantages and disadvantages of both

methods from a theoretical as well as a numerical point of view.

Advantages of the Gauss-Galerkin method

 

1. When the truncated system of equations for the nodes and the weights admits
a solution on a t interval of the form [0,T], where T 6 (0,00), then the

approximate solution converges to the law of X, in distribution.

2. The approximation in part 1 allows any functional of the process X to be

approximated.

3. Using the numerical algorithm we proposed in Chapter 4, the Gauss-Galerkin
method proved to be numerically efficient.

Disadvantages of the Gauss-Galerkin method

 

1. The Gauss-Galerkin method fails to work when the initial distribution is

degenerate.

2. The approximate solution for the process moments has no known rate of

convergence.

3. The existence of a unique solution for the truncated system of equations for
the nodes and weights (equivalently for the truncated system of moments
when the drift and the square of the diffusion are polynomials) is only true
on a t interval of the form [0,T(n)] where T(.) is a function of n. There is
no available theoretical proof for a global solution, i.e. T(n) E T, Vn E N.

However, for some models, this can be established numerically.

Advantages of the power series method

 

45

1. The truncated system of equations for the moments has a unique solution on
any fixed t interval.

2. The approximate solution converges to moments of X uniformly on any in-

terval [0,T], for all T 6 (0, oo).

3. The expression of the approximate solution can be explicitly found as a func-

tion of t, for all t 6 (0,00).

Disadvantages of the power series method

 

1. It is not clear that the truncated system of equations for the moments yields
functions which are moments of a probability measure. This does not make

approximations of functionals of the process possible.

2. As in the Gauss-Galerkin method, there are no available rates of convergence

to the steady state solution of the stochastic differential equations.

4.3 Steady State Analysis

In this section we study the behavior of the process X as t —» 00. According to

Theorem 3.1 (Chapter 3), we have VT > 0, and We 6 N,
lim sup IMk(t) - mk(n,t)| = 0.
""°° 05:57"
For any given k E N we can uniformly approximate the k-th moment by the

sequence {mk(n, .), n 2 1} on [0,T]. When we let t increase to infinity, we have
mk(n,t) —r 00.

Therefore to analyze the steady state solution, a single polynomial of degree n
(no matter how large n is) will not provide useful information when t —; 00. To
illustrate this we choose the initial distribution to be uniform and fix n = 7. Using

the power series method described earlier, the following results are derived for
Model 3,

7
m2(7, t) = 2 c,-t',
i=0

46

the coefficients c,, for i = 0, 1, ..., 7 are as above
The rest of the functions m3(7, t), m4(7, t), ..., m9(7, t) are easily obtained from

the function m2(7, t). Moreover, it is clear that
111.12, m2(7,t) = 00.

Based on the power series method, the following computational algorithm yields
useful information about the steady state solution of the process X. The results are
summarized in the appendix. We now give an outline of this algorithm designed

for the study of the steady state solution.

Computational Algorithm

 

1. Given any T 6 (0,00), let h = %, where N E N.

2. Divide the interval [0, T] into N sub-intervals of length h, to get the partition
[O,h1)[h1,h2), . . ., [hN-1,hN], where h,- =1h, 0 S13 N — l.

3. Fix a low degree p (say p = 1 or p = 2).

4. Define the following functions
- p (i)
m(k, 1,t) = z cj (k)(t — h,-_1),
i=0

where t E [h,-_1,h,-) and 1 S i _<_ N — 1, and m(k,i,t) will represent the
piecewise linear (if p = 1) or quadratic (if p = 2) approximation of the k-th
moment in the i-th sub-interval. For p = l, the functions {m(k, i,t),0 5 i g

N — 1} are computed as follows.
5. For the sub—interval [0, h), we have
m(2.1,t) = of.” + c‘.”t.

where c9) = M2(0) and c9) = Mi’,(0). This yields an approximation for M2(t),
t E [0, h). Note that we used two initial conditions. In general, to calculate
m(k, l, t), for each k, 2 initial conditions are needed. Computing 0 functions

will require oz+1 initial conditions.

47

6. For the sub-interval [h, 2h), we have
m(2, 2,.) = ct” + 4% — h),

where 682) = m(2,1,h) and C(12) = 02(m(2,1,h) —- m(3, 1, h)). Note that com-
puting more and more functions will involve more and more initial conditions.
As we march on from sub-interval to sub-interval, more initial conditions are
needed.

7. To study the steady state solution for M2(t) for example, we observe that
(lim CI = 0,
this yields the steady state value (for M2)

lim c0.

l—boo

A similar behavior is observed for higher moments.

Recall that we considered the operator L over the interval [0,1] with its domain
containing all the monomials. This assumption holds for all three models. The
boundary 0 and 1 may be classified as ”regular”, ”exit”, ”entrance” or ” natural”
(see Appendix), depending on the behavior of the functions a and b involved in the
stochastic differential equation. Using the boundary classification, we can easily

show the following:

In Model 1, both 0 and l are ”natural” boundaries. In Model 2, both 0 and 1
are ”exit” boundaries. In model 3, 0 is a ”natural” boundary and 1 is an ”exit”
boundary. In model 4.1, both 0 and l are ”entrance” boundaries. In model 4.2,
both 0 and 1 are ”exit” boundaries. According to Ethier and Kurtz, the appropriate

domain of the operator L is given by
D(L) .—. {<15 6 C[0,1] 0 02(0, 1) : L(qb) e C[0,1]},
at an entrance boundary or at a natural boundary, and

D(L) = {as e 010.11n 010,1) = Les) 6 C[O. 11.12:; Lab) = 0}.

48

at an exit boundary r. Moreover, all of our assumptions (assumptions A1 through
A5) are satisfied by these models.

In Tables 1 through 3, numerical results for the moments mo, m1, . . . , m5 are pre-
sented for several t values until the steady state solution is reached. Table 1 is for
model 1, Table 2 is for model 2, and Table 3 is for model 3. The initial distribu-
tion considered is the uniform on the interval (0,1). Since the function a vanishes
identically on the state space [0, 1] in all three models, it follows using system (3.9)
that both mo and m, are constants in t, and according to the initial distribution

(uniform) their values are 1 and 0.5, respectively.

In Model 1, both boundaries are natural and the law of X, is absolutely contin-
uous with respect to the Lebesgue measure on [0,1]. In this model the density
p(oo, y) E 0, Vy 6 (0,1), is clearly a steady state solution with the density escaping
to the boundary as t becomes arbitrary large. Our computations show that the
steady state solution is reached when t is approximately 50. The results in Table 1
show that the steady state values for the moments are approximately equal to 0.5.
This implies that the steady state distribution measure for the process X is the
measure $760 + 61). This is in agreement with the exact solution of the model (see
Appendix). Table 1.1 describes the approximate values of the first few moments

through the power series method until steady state is reached.

In Model 2, both boundaries are exit and the law of X, is not absolutely con-
tinuous with respect to the Lebesgue measure on [0,1]. However, the absolutely
continuous part of the law of X, has support (0, 1). The singular part of the law
of X, has a point mass at 0 and a point mass at 1. It is important to note that
the Radon-Nikodym derivative p(t, y) of the absolutely continuous part satisfies the
Fokker-Planck equation. In this model, for each y 6 (0,1), the function p(t, y) —i 0,
as t —» 00. The steady state solution for this model can be deduced from the be-
havior of the moments as t -> 00. Our computations show that the steady state
solution is reached when t is approximately 12. The results in Table 2 show that

the steady state values for the moments are approximately equal to 0.5. This im-

49

plies that the steady state distribution measure for the process X is the measure
-;-(60 + 61). This is in agreement with the exact solution. Table 2.1 describes the
approximate values of the first few moments through the power series method until

steady state is reached.

In Model 3, 0 is a natural boundary and 1 is an exit boundary. In this model the
law of X, is again not absolutely continuous with respect to the Lebesgue measure
on [0,1]. The singular part of the law of X, has a point mass at 1. It is important to
note that the Radon-Nikodym derivative p(t, y) of the absolutely continuous part
satisfies the Fokker-Planck equation. In this model the function p(t,y) -» 0, as
t -+ 00. The steady state solution for this model can be deduced from the behavior
of the moments as t -> 00. The results in Table 3 show that the steady state values
for the moments are approximately equal to 0.5. This suggests that the steady
state distribution measure for the process X is the measure lz-(bo + 61). Our compu-
tations show that the steady state solution is reached when t is approximately 20.
The steady state analysis for this model cannot be checked against an exact solu-
tion since no exact solution is known. Table 3.1 describes the approximate values of

the first few moments through the power series method until steady state is reached.

In Model 4.1, both 0 and 1 are entrance boundaries. In this model the law of
X, is absolutely continuous with respect to the Lebesgue measure on [0,1]. The
Radon-Nikodym derivative p(t, y) of the law of X, satisfies the Fokker-Planck equa-
tion. The steady state solution for this model can be found by solving the equa-
tion L*p(oo, y) = 0. Using both the Gauss-Galerkin method and the power series
method, we can recover the steady state moments. By applying the power series
method to this model we observe that the steady state solution is reached at ap-
proximately t = 6. By applying the Gauss-Galerkin method, we concluded that the
steady state solution is achieved at approximately t = 5. The exact steady state
moments are very close to both the approximate steady state moments obtained
through the power series method and to the approximate steady state moments ob-

tained through the Gauss-Galerkin method. Table 4.1 describes the approximate

50

values of the first few moments through the power series method until steady state
is reached. Table 4.1.1 is a summary of the true steady state moments and Table

4.1.2 shows the approximate steady state moments.

In Model 4.2, both 0 and l are exit boundaries. In this model the law of X, is not
absolutely continuous with respect to the Lebesgue measure on [0,1]. The singular
part of the law of X, has a point mass at 0 and 1. It is important to note that
the Radon-Nikodym derivative p(t, y) of the absolutely continuous part satisfies the
Fokker-Planck equation. In this model the function p(t,y) —> 0, as t —> 00. The
results in Table 3 show the approximate values of the moments for several values
of t. The exact solution is not known. Therefore we did not compare our results
with the values of the true moments. Table 4.2.1 is a summary of the values of the
first few moments obtained through the power series method. Table g.4.2 shows
the approximate values of the first few moments through the power series method.
The analysis of the steady state behavior of this process is work in progress and

will be included in a later publication.

4.4 Discussion

This section provides concluding remarks and an outline of some suggestions for
future work. One open problem is the approximation of joint distributions of the
process X. A more difficult task is the approximation of the transition probabil-
ity function. An examination of the conditional moments may lead to some useful
ideas for the approximation of the transition density function. We also outline some
problems which are still unresolved for both the Gauss-Galerkin method and the
power series method. An important task is the development of a numerical method
which provides an approximation of the transition probability function and there-
fore leading to an approximation of the joint distribution of the diffusion (since any

joint distribution of the process X is in terms of the transition probability function).

We propose an idea that might lead to some new development in this direction. We

51

note that neither the power series method nor the Gauss-Galerkin method offers

this kind of an approximation. An outline of our idea is as follows.

The stochastic differential equation 3.1 has a unique solution X, which is a time
homogeneous Markov process. Moreover, assume the density p(t, 2:, y) of the tran-

sition probability function satisfies the following initial boundary value problem

Eat—5,5439- = L.(p<t.4.y>). (4.10)
p(0,w.y) = 6(x—y), (4.11)
p(th) = ¢(t.y). (4.12)
p(tmy) = 1120.21). (4.13)

where the state space of the process X is one of the following (l,r), [l,r], (l,r],
or [1, r). The functions a) and 1b are given. Because of the lack of smoothness of
the initial function p(0, r, y), the above initial boundary value problem may be an
ill-posed problem and hence difficult to handle.

A somewhat easier problem for the conditional moments of the process X can be

derived as follows.

The i - th conditional moment of the process is defined by

0.0.30) = [,y‘p(t,x,y)dy,w 6 W U 0.

(assuming the integrals exist). In Equation 4.10 we multiply both sides of the

equation by y‘, then integrate with respect to the variable y.

Then the following initial boundary problem can be derived

0C5“, 27)
T L,(C,-(t, r)), (4.14)

C,‘(0,£L‘) = (1}., (4.15)

52

C,(t, I)

[1' W. y)y‘dy. (4-16)
C.(t. r) = [Ir W. y)y‘dy- (4-17)

Note that this system does not suffer the lack of smoothness of the initial condi-
tion. Moreover, the system for the conditional moments may be easier to solve.
Moreover, the process moments may be obtained by integration of the conditional
moments. A much more important task is the approximation of the distribution of

X, be a sequence of discrete probability measures for each value of t.

The knowledge of the first it conditional moments may lead to the construction of
probability measures which we conjecture will converge to the transition probability
measure. More work along the above ideas is in progress. Another problem which
more complicated is the two dimensional problem. In some areas of engineering
sciences, stochastic differential equations where the solution is a two dimensional
stochastic process present some challenging questions. One of the most difficult
problems is the approximation of the two dimensional distributions involved. More

work needs to be done in this direction.

APPENDIX

53

APPENDIX

In this appendix we present the following items. Section A.1 contains a summary
of the figures mentioned in Chapter 4. Section A.2 contains a summary of the
tables mentioned in Chapter 4. Section A.3 contains a discussion on boundary

classifications.

A.1: Figures 1 and 2 are results of an application of the moment neglect method
to Model 3. Figure 1 shows a plot of the initial moments M2(0), M3(0), . . . , M3(0)
against the order of the moments. Figure 2 shows a plot of the approximate mo-
ments (according to the moment neglect method) against the order of the moments,
for t = 0.001. It can be seen from Figure 2 that the results are misleading since the
approximate value of M6(0.001) clearly exceeds 1. This contradicts the fact that
the state space of the diffusion under consideration is [0, 1]. Figures 1.1, 1.2, 1.3,
1.4, 1.5, 1.6, and 1.7 represent the exact transition probability density for model 1
plotted against the forward variable y, for several values of t until the exact steady
state is reached. Figures 1.2 through 1.4 show that for t between 0 and 10 the
density is unimodal. For t = 20, approximately, the density is shown in Figure 1.5
to be bimodal. The density eventually accumulates on the boundaries 0 and 1 at
approximately t = 50. Figure 1.7 shows that both accumulations at 0 and 1 are
approximately equal to 0.5. This is in agreement with both the Gauss-Galerkin
method and the power series method. Figures 2.1, 2.2, 2.3, 2.4, and 2.5 represent
the exact Radon-Nikodym derivative of the absolutely continuous part of the tran-
sition probability measure with respect to the Lebesgue measure. This function is
plotted against the forward variable y, for several values of t until the steady state
is reached. Figure 2.5 shows for t = 12, approximately, the steady state value of
the Radon-Nikodym derivative of the transition completely leaks out through both
boundaries of the state space [0,1]. This ”numerical leakage” is due to the fact that

both boundaries 0 and 1 are exit boundaries.

A.2: Table g.2 represents the approximate moments for Model 2 according to the

Gauss-Galerkin method. In Table g.2, m1,m2, . . . ,m5 represent the approximate

54

values of the exact moments, M1,M2, . . . ,M5, of the process. The approximate
steady state moments stabilize around the value 0.5 suggesting that the steady state
density has a point mass at 0 and a point mass at 1, both equal to 0.5. Table g.3
represents the approximate moments for Model 3 according to the Gauss-Galerkin
method. A similar behavior to the one in Table g.2 is observed. Since the exact
moments are not known for this model, we cannot make quantitative comments on
the accuracy of the results. However, the power series method is in agreement with
these results (see Table 3.1). Table g.4.1, g.4.2 represent the approximate moments
for Model 4.1 and Model 4.2, respectively, according to the Gauss-Galerkin method.
Both tables are in agreement with the corresponding results obtained by using the
power series method (see Table 4.1 and 4.2.1). Table 1.1, 2.1, 3.1, 4.1, and 4.2.1
represent the power series approximate moments for Models 1, 2, 3, 4.1, and 4.2,
respectively. In these tables, the values of the functions m1, m2, . . . , m5 are given
for several values of t. The functions m1,m2, . . . ,m5 represent the approximate
values of the exact moments, M1, M2, . . . , M5, of the process under consideration.
From Table 1.1 it can be seen that the steady state solution for Model 1 is reached
for t = 55, approximately. We also note that the approximate values of the mo-
ments stabilize around 0.5 suggesting that the steady state density has has a point
mass at 0 and a point mass at 1, both equal to 0.5. This is in agreement with the
Gauss-Galerkin method. A similar behavior can be observed through Tables 2.1,
and 3.1 for models 2, and 3, respectively. Table 4.1 shows that the approximate
steady state solution is reached at t = 20, approximately. This agrees with the
exact steady state solution. Tables 4.1.1, and 4.1.2 show this agreement. Table
4.2.1 is for Model 4.2. It shows the power series approximate values for several
values of t. For this model we do not have a steady state analysis at this point.
We have compared our results (both through the Gauss-Galerkin method and the
power series method) with the results obtained by J. Keller and R. Voronka and

we observed a close agreement.

55

A.3. Boundary Classification
Let X be a diffusion process with state space [0,1]. According to Feller [15], the

boundaries 0 and 1 can be classified as follows.

1. Regular.

2. Exit.

3. Entrance.

4. Natural (or Feller).

We shall present the formal definitions of the different types of boundaries. These
definitions can be found in Karlin and Taylor, and Feller. The scale function is
defined by

3

0(4) = j, q(y)dy,

0

where
(M) = epr- A: 20(n)b‘2(n)dnl.Vy 6 (0.1).

The numbers 1:0 and yo are arbitrary, but fixed. Their particular choice is of no

relevance to the analysis below. The scale measure is defined by
u(d:r) = q(a:)da:.

The speed measure is defined by
duh?) = m($)dx,

where m(r) = W'

Using the scale and the speed measures, the following important quantities are

defined. For the left boundary point 0, we have

2(0) = [4((4,4))q(e)44,

and
9(0) = f,” 441:, 41)m(c)44.

The following summarizes the boundary classifications.

56
0 is regular if 2(0) < 00 and 9(0) < oo.
0 is an exit if 2(0) < 00 and 9(0) = 00.
0 is an entrance if 2(0) = 00 and 9(0) < oo.
0 is natural if 2(0) = 00 and 9(0) < 00.

The boundary behavior at 1 is defined similarly.

57

Moment Neglect Method (Model 3)

 

Figure 1: Moment Neglect, Time=0

 

 

 

 

10 Figure 2: Moment Neglect, Time=0.001

 

N
U10U'1U'IU'I

-7.S'

 

A A A A A A
-

 

 

Model 1 (exact solution)

58

 

Figure 1.1, Model 1

 

O
O
N
(11

O
O
..s
U1

P(0.1.0.S.y)

0.005

O
O
N

0.01)

 

D
>#

 

M

v

 

1
L 1
A A

 

0

0.2

:4

o:

6

0

Figure 1.2, Model 1

.8

1

 

 

ONIhO‘mt-‘Nlb

f

 

 

 

 

59

Figure 1.3. Model 1

 

 

 

P(3.0 5.y)
c: c: c: c: c: c:
OHlvlaJcbuim

 

 

Figure 1.4, Model 1

 

T v f 1

0.145
0.12»

.y)
0
H

-o.08»
‘0.06'

P(10 0 S

O
O
A

0.02-

 

 

 

 

Figure 1.5, Model 1

 

v r r v V

 

 

P(20,0.5.Yl
<3 :1 <3 c: <3
w

 

 

 

60

Figure 1.6. Model 1

 

N

...:
OU'Il-‘UINUIW

P(35,0.5.yl

O

 

V V v

 

Or

0:2 0:4 0:6 0:

Figure 1.7, Model 1

 

h) lb
0 O

N
O

 

P(S0.0.S.y)

..s
O

 

a v M Y

 

 

0:2 0:4 0:6 ole

 

 

 

61

Model 2 (exact solution)

Figure 2.1, Model 2, Exact Solution

 

f

i

.S.y)
O

P(l,

 

 

 

 

Figure 2.2, Model 2, Exact Solution

17

 

T f 1

0.2)

0.15»

.5,y)

0.1)

P(4,

0.05'

 

 

 

62

Figure 2.3, Model 2, Exact Solution

0.005)

 

 

r v v v *v

 

 

 

Figure 2.4, Model 2, Exact Solution

0.01v

0.008)
50.006»
30.004.

0.002)

 

 

——v v w

 

 

0:2 0:4 0:6 0:8 1

Figure 2.5, Model 2, Exact Solution

 

0.4)

0.2)

.5,y)

 

P(12.

 

 

 

0:2 0:4 0:6 ole i

63

Gauss Galerkin Method, Model 2

 

Approx. Moments, Model 2, Table g.2

 

t

m1

m2

m3

m4

m5

 

 

0.000000
2.000000
8.000000
30.00000

 

0.500000
0.497900
0.497890
0.497690

 

0.333333
0.475130
0.497830
0.497680

 

0.250000
0.467500
0.497800
0.497680

 

0.250000
0.456920
0.497780
0.497680

 

0.166667
0.452350
0.497770
0.497680

 

Gauss Galerkin Method, Model 3

 

 

Approx.

Moments, Model 3, Table g.3

 

t

m1

m2

m3

m4

m5

 

 

0.000000
2.000000
8.000000
30.00000

 

0.500000
0.501030

0.501000.

0.500840

 

0.333333
0.427030
0.484760
0.498700

 

0.250000
0.401640
0.481690
0.498600

 

0.250000
0.389110
0.480380
0.498560

 

0.166667
0.381730
0.479660
0.498540

 

 

 

64

Gauss Galerkin Method, Model 4.1

 

 

Approx. Moments, Model 4.1, Table g.4.1

 

t

m,

m2

m3

m4

m5

 

 

0.000000
2.000000
4.000000
6.000000
10.00000
20.00000

 

0.500000
0.501100
0.490003
0.498879
0.498888
0.502010

 

0.333333
0.290225
0.299002
0.301100
0.301230
0.300400

 

0.250000
0.200337
0.200022
0.202300
0.199050
0.203201

 

0.167000
0.143231
0.142860
0.142863
0.142863
0.142863

 

0.200000
0.107517
0.107139
0.107300
0.107566
0.107102

 

Gauss Galerkin Method, Model 4.2

 

 

Approx. Moments, Model 4.2, Table g.4.2

 

t

m1

m2

m3

m4

m5

 

 

0.000000
1 .000000
2.000000
3.000000
4.000000
5.000000
6.000000

 

0.500000
0.520564
0.599000
0.616667
0.645500
0.670098
0.703976

 

0.333333
0.377767
0.412280
0.442779
0.478420
0.511572
0.541934

 

0.250000
0.279300
0.319802
0.344444
0.379263
0.412804
0.444842

 

0.20000
0.225088
0.251222
0.284367
0.314702
0.348146
0.378972

 

0.166667
0.189048
0.214008
0.241379
0.270711
0.305273
0.335410

 

 

 

65

Power Series Method, Model 1

 

 

Approx. Moments, Model 1, Table 1.1

 

t

m1

m2

m3

m4

m5

 

 

0.000000
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
0.700000
0.800000
0.900000
1.000000
5.000000
12.00000
35.00000
55.00000

 

0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000

 

0.333333
0.336667
0.339905
0.343052
0.346114
0.349093
0.351994
0.354819
0.357573
0.360257
0.362875
0.394726
0.474390
0.492492
0.495014

 

0.250000
0.255000
0.259857
0.264579
0.269171
0.273640
0.277991
0.282229
0.286359
0.290386
0.294313
0.353642
0.456862
0.489186
0.494468

 

0.250000
0.205714
0.21 1286
0.216719
0.222021
0.227194
0.232243
0.237174
0.241989
0.246693
0.251289
0.325734
0.443768
0.486220
0.492990

 

0.166667
0.172619
0.178452
0.184168
0.189767
0.195251
0.200623
0.205886
0.211040
0.216089
0.221035
0.304407
0.432859
0.483435
0.491564

 

 

66

Power Series Method, Model 2

 

 

Approx.

Moments, Model 2, Table 2.1

 

t

m1

m2

m3

m4

m5

 

 

0.000000
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
0.700000
0.800000
0.900000
1.000000
2.000000
3.000000
4.000000
5.000000
6.000000
7.000000
8.000000
10.00000
12.00000

 

0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0,500000
0.500000
5.000000
5.000000
5.000000
5.00000

5.000000
5.000000
5.000000

 

0.333333
0.349194
0.363545
0.376530
0.388280
0.398912
0.408531
0.417236
0.425112
0.432238
0.438687
0.477444
0.491702
0.496947
0.498877
0.499587
0.4998480
0.499944
0.499973
0.499838

 

0.250000
0.273791
0.295317
0.314795
0.332420
0.348367
0.362797
0.375854
0.387668
0.398358
0.408030
0.466166
0.487553
0.495421
0.498316
0.499380
0.499772
0.499916
0.499960
0.499757

 

0.200000
0.228549
0.254381
0.277755
0.298904
0.318041
0.335357
0.351024
0.365201
0.378029
0.389636
0.459399
0.485064
0.494505
0.497979
0.499256
0.499726
0.499899
0.499952
0.499709

 

0.166667
0.198388
0.227090
0.253061
0.276560
0.297823
0.317063
0.334472
0.350224
0.364477
0.377374
0.454888
0.483404
0.493895
0.497754
0.499174
0.499696
0.499888
0.499947
0.499677

 

 

67

Power Series Method, Model 3

 

 

Approximate Moments, Model 3, Table 3.1

 

t

m1

m2

m3

m4

m5

 

 

0.000000
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
0.700000
0.800000
0.900000
1 .000000
3.000000
12.00000
20.00000

 

0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000

 

0.333333
0.341667
0.349333
0.356417
0.362985
0.369094
0.374793
0.380122
0.380122
0.389807
0.394220
0.420367
0.482793
0.499200

 

0.250000
0.265000
0.278500
0.290736
0.301889
0.312106
0.321503
0.330179
0.330179
0.345678
0.352630
0.404972
0.480233
0.498951

 

0.167000
0.220000
0.237714
0.253557
0.267834
0.280781
0.292584
0.303394
0.303394
0.322506
0.330997
0.382445
0.480088
0.498515

 

0.200000
0.190476
0.211310,
0.229762
0.246256
0.261109
0.274568
0.286828
0.286828
0.308355
0.317859
0.384063
0.479000
0.498095

 

 

68

Power Series Method, Model 4.1

 

 

Approx. Moments, Model 4.1, Table 4.1, u=v=h=0.5, and s=0

 

t

m1

m2

m3

m4

m5

 

 

0.000000
0.100000
0.300000
0.500000
0.700000
0.900000
1.100000
1.300000
1.700000
1.900000
2.000000
4.000000
6.000000
8.000000
10.00000
12.00000
14.00000
16.00000
18.00000
20.00000

 

0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000

 

0.333333
0.325960
0.315746
0.309550
0.305792

0.303513

0.302131
0.301292
0.300475
0.300288
0.300225
0.300002
0.300000
0.300000
0.300000
0.300000
0.300000
0.300000
0.300000
0.300000

 

0.250000
0.238940
0.223618
0.214325
0.208689
0.205270
0.203196
0.201939
0.200713
0.200433
0.200337
0.200002
0.200000
0.200000
0.200000
0.200000
0.200000
0.200000
0.200000
0.200000

 

0.167000
0.186912
0.169294
0.158822
0.152523
0.148716
0.146409
0.145011
0.143650
0.143338
0.143231
0.142860
0.142863
0.142863
0.142863
0.142863
0.142863
0.142863
0.142863
0.142863

 

0.200000
0.152380
0.133871
0.123180
0.116827
0.113006
0.110696
0.109297
0.107935
0.107624
0.107517
0.107145
0.107142
0.107142
0.107142
0.107142
0.107142
0.107142
0.107142
0.107142

 

 

69

Steady State, Model 4.1

 

Steady State, Model 4.1, Figure 4.1

 

6y(1-y)

O O O O
o 0 O O o O
O N ab 0‘ (D H N 1b

 

#1

v7

 

 

 

 

True Steady State Moments, Model 4.1, Table 4.1.1

 

m,

7H2

1123

m4

m5

 

 

0.500000

 

0.300000

 

 

0.200000

0.142857

 

0.107143

 

 

Approx. Steady State, Model 4.1, Table 4.1.2, t=20

 

772,

m2

m3

m4

m5

 

 

0.500000

 

0.300000

 

0.200000

 

0.142863

 

0.107142

 

 

 

70

Power Series Method, Model 4.2

 

 

Approximate Moments, Model 4.2, Table 4.2.1

 

t

m,

m2

m3

m4

.m5

 

 

0.000000
1.000000
2.000000
3.000000
4.000000
5.000000
6.000000

 

0.500000
0.540000
0.578000
0.612767
0.643648
0.670598
0.693976

 

0.333333
0.368867
0.405886
0.442871
0.478420
0.511572
0.541934

 

0.250000
0.279300
0.311312
0.345041
0.379263
0.412804
0.444842

 

0.20000
0.225057
0.252947
0.283167
0.314902
0.347146
0.378972

 

0.166667
0.189048
0.214008
0.241379
0.270711
0.301273
0.332219

 

 

BIBLIOGRAPHY

71

BIBLIOGRAPHY

Arnold L. (1974). Stochastic Differential Equations, Theory and Applications. Wi-
ley Interscience.

Dawson, D. A. (1980). Galerkin Approximation of Nonlinear Markov Processes,
Statistics and Related Topics. Proceedings of the International Symposium on
Statistics and Related Topics held in Ottawa, Canada, May 5—7, 1980. North Hol-
land.

Ethier, Stewart N. and Kurtz, Thomas G., (1986). Markov Processes, Characteri-
zation and Convergence. John Wiley and Sons, New York.

Ewens, Warren J. (1979). Mathematical Population Genetics. Biomathematics,
Volume 9, Springer-Verlag, New York.

Feller, W. (1952). The Parabolic Differential Equation and the Associated Semi-
Group of Transformation, Ann. Math., Vol. 55, pp. 468-519.

HajJafar, A. (1984). On the Convergence of the Gauss Galerkin Method for the
Density of Some Markov Processes. Ph.D. Thesis, Department of Mathematics,
Michigan State University.

HajJ afar, A., Salehi, H., and Yen, D. (1992). On a Class of Gauss-Galerkin Methods
for Markov Processes, Proceedings of the Work Shop on Nonstationary Stochastic
Processes and Their Applications, Edited by Miamee, A. G., Hampton, Virginia,
August 1-2, 1991, pp. 126—146.

Karlin, S. and Taylor, H. M. (1981). A Second Course in Stochastic Processes.
Academic Press.

Keller, J. and Voronka, R. (1975). Asymptotic Analysis of Stochastic Models in
Population Genetics. Mathematical Biostatistics, Vol. 25, pp. 331-362.

Kimura, M. and Crow, J. (1970). An Introduction to Population Genetics. Harper
and Row, New York.

Krylov, V. I. (1962). Approximation Calculation of Integrals. Macmillan, New
York. (Translated from the first Russian edition (1959) by Stroud, A. H.)

Shohat, J. A. and J. D. Tamarkin (1943). The Problem of Moments, A.M.S., New

72

York.

Smirnov, V. I. (1964). A Course of Higher Mathematics. Pergamon Press, Oxford,
New York.

Stoer, J. and Bulirsch, R. (1980). Introduction to Numerical Analysis. Springer-
Verlag, New York.

Stroud, A. H. ( 1974). Numerical Quadrature and Solution of Ordinary Differential
Equations. Springer-Verlag, New York.

Szarski, J. (1965). Differential Inequalities. Warszawa.

Wentzell, A. D. (1981). A Course in the Theory of Stochastic Processes. McGraw
Hill, New York.

111011111111)“