\magnification = \magstephalf
\hsize=14truecm
\hoffset=1truecm
\parskip=5pt
\nopagenumbers
\font\eightrm=cmr8 \font\eighti=cmti8 \font\eightbf=cmbx8
\headline={\ifnum\pageno=1 \hfill\else%
{\tenrm\ifodd\pageno\rightheadline \else
\leftheadline\fi}\fi}
\def\rightheadline{EJDE--1995/15\hfil Approximate General Solution 
\hfil\folio}
\def\leftheadline{\folio\hfil Kazuo Amano 
\hfil EJDE--1995/15}
\voffset=2\baselineskip
\vbox {\eightrm\noindent\baselineskip 9pt %
 Electronic Journal of Differential Equations,
Vol. {\eightbf 1995}(1995) No.\ 15, pp. 1--14.\hfill\break
ISSN 1072-6691. URL: http://ejde.math.swt.edu (147.26.103.110)\hfil\break
telnet (login: ejde), ftp, and gopher access:
 ejde.math.swt.edu or ejde.math.unt.edu}
\footnote{}{\vbox{\hsize=10cm\eightrm\noindent\baselineskip 9pt %
1991 {\eighti Subject Classification:} 35K65, 65M99.
\hfil\break
{\eighti Key words and phrases:} Degenerate parabolic, numerical-symbolic
method.
\hfil\break
\copyright 1995 Southwest Texas State University  and
University of North Texas.\hfil\break
Submitted February 20, 1995. Published October 20, 1995.} }

\bigskip\bigskip

\centerline{APPROXIMATE GENERAL SOLUTION OF DEGENERATE}
\centerline{PARABOLIC EQUATIONS RELATED TO POPULATION GENETICS}
\smallskip
\centerline{Kazuo Amano}
\bigskip\bigskip

{\eightrm\baselineskip=10pt \narrower
\centerline{\eightbf Abstract}
The author constructs an approximate general solution to a
degenerate parabolic equation related to population genetics
and implements a computational procedure. 
The result gives a theoretical foundation to the computer algebraic 
approach for degenerate partial differential equations and introduces
a new numerical symbolic hybrid method.
The techniques are likely to have wide applicability,
since the key idea of the algorithm is a rearrangement of 
the finite difference method.
\bigskip}

\def\Proof{\smallbreak\noindent{\bf Proof.\enspace}}
\def\QED{\ \hbox{\vrule height 1.6ex width 1ex depth .3ex}\ }
\def\pd{\partial}
\def\ds{\displaystyle}

\bigbreak
\centerline{\bf \S 1. Introduction} \medskip\nobreak

As is well-known, if a given partial differential equation is very simple,
we can compute its general solution with arbitrary functions
by using arithmetic and elementary calculus.
However, most equations require
hard and abstract mathematical technicalities.
An explicit and concrete representation of the solution
may turn out to be utterly beyond our reach.
It seems that
researchers have already given up constructing explicit general solutions;
they are either trying to find solutions in abstract function spaces
or working out numerical algorithms.

In this paper we shall show
new possibilities for {\it approximate general solutions\/};
though an explicit representation is already at deadlock,
an approximate one is able to break through obstacles
and gives a new viewpoint.
In fact, we prove that,
for a certain initial value problem,
there exists a simple algebraic representation
of an approximate general solution,
i.\ e., a symbolic combination of additions, subtractions and multiplications
of initial data solves the problem.
Our procedure of construction of a general solution
is quite different from classical ones;
we use a new type of {\it numerical-symbolic hybrid method\/}.
Our numerical-symbolic hybrid computation
totally depends on LISP
and its result is expressed in C language,
since the size of desired formula is more than 5.4M bytes.
Such a formula is too big for classical pen and paper calculation.
It is to be noted that a remarkably fast algorithm is derived
from our formula of approximate general solution.

The purpose of this paper is to construct an approximate general solution
of the initial value problem for the degenerate parabolic equation
$$
{\pd u\over\pd t}
	={1\over2}{\pd^2\over\pd x^2}\bigl(V(x)u\bigr)
	-{\pd\over\pd x}\bigl(M(x)u\bigr)\qquad (t>0,\ 0<x<1)\ ,
\leqno(1.1)
$$
which appears in population genetics theory.
$\,\ds V(x)={x(1-x)\over2N}\,$,
$\,M(x)\,$ is a polynomial of $\,x\,$
, and $\,N\in\{1,2,3,\cdots\}\,$ denotes the given population number.
For a certain class of coefficients $\,M(x)\,$,
Crow and Kimura obtained several explicit representations of solutions of (1.1)
by using special functions (cf.\ [2]).
However, their technicality is too artistic and too delicate to generalize;
in fact, if you perturb the function $\,M(x)\,$ slightly,
their method simply breaks down.
In this paper, 
we shall show that we can construct an approximate general solution of (1.1)
for any given $\,M(x)\,$.
Our techniques are applicable to
the broader class of differential equations for which
finite difference method works,
since the essential part of our method is nothing more than
a rearrangement of a finite difference scheme (Lemma 2.3).

Section 2 is devoted to preliminary lemmas.
In Section 3, we construct an approximate general solution of (1.1)
and translate its procedure into symbolic list operations,
which is expressed in PASCAL-like REDUCE language ([4]).
We make some numerical experiments in Section 4.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%					%
%		Section 2		%
%					%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigbreak
\centerline{\bf \S 2. Preliminaries} \medskip\nobreak

We simplify the partial differential equation (1.1)
by performing a certain transformation (Lemma 2.2)
and give an approximate representation formula (Lemma 2.3).

\proclaim {Lemma 2.1}.
For any nonnegative integer $\,n\,$
and for any $\,C^{n+1}\,$ smooth function $\,u(t,x)\,$,
we have
$$
\eqalign{
u(t+h,x+k)
&=\sum_{\nu =0}^n{1\over\nu !}
	\Bigl(h{\pd\over\pd t}+k{\pd\over\pd x}\Bigr)^{\nu }
	u(t,x)\cr
&+\int_0^1{(1-\theta )^n\over n!}
	\Bigl(h{\pd\over\pd t}+k{\pd\over\pd x}\Bigr)^{n+1}
	u(t+\theta  h,x+\theta  k)\,d\theta \cr
}
\leqno(2.1) $$

\Proof
We fix $\,(t,x)\,$ arbitrarily and put
$$
F(\theta )=u(t+\theta  h,x+\theta  k)\qquad (0\leq \theta \leq 1)\ .
$$
We note that
$$
u(t+h,x+k)=F(1)=\sum_{\nu =0}^n{1\over\nu !}\,F^{(\nu )}(0)
	+\int_0^1{(1-\theta )^n\over n!}\,F^{(n+1)}(\theta )\,d\theta \ .
\leqno(2.2) $$
When $n=0$, (2.2) is nothing more than fundamental theorem of calculus.
Integration by parts gives
$$
\eqalign{
&\int_0^1{\ \ (1-\theta )^{n-1}\over(n-1)!}
	\,F^{(n)}(\theta )\,d\theta \cr
&=\int_0^1\Bigl\{-{(1-\theta )^n\over n!}\Bigr\}'
	\,F^{(n)}(\theta )\,d\theta \cr
&={1\over n!}\,F^{(n)}(0)
	+\int_0^1{(1-\theta )^n\over n!}
	\,F^{(n+1)}(\theta )\,d\theta \ .\cr
}
\leqno(2.3)$$
Hence, (2.1) is proved.
\QED

Let us define a differential operator $\,A\,$ by
$$
Au={1\over2}{\pd^2\over\pd x^2}\bigl(V(x)u\bigr)
	-{\pd\over\pd x}\bigl(M(x)u\bigr)\ ,
\leqno(2.4) $$
where
$\ds\,V(x)={x(1-x)\over2N}\,$,
$\,M(x)\,$ is a polynomial of $\,x\,$ with real coefficients,
and $\,N\equiv N_e\in\{1,2,3,\cdots\}\,$ denotes the effective 
population number. 
As is well-known, (2.4) plays an essential role in
population genetics theory (cf.\/, for example, [2]).
In fact, if we assume that
a pair of alleles $\alpha _1$ and $\alpha _2$ are segregating
in a Mendelian population,
then the probability density $u(t,x)$ of frequency of $\alpha _1$
at the $\,t\,$th generation satisfies the partial differential equation
$$
{\pd u\over\pd t}=Au\hbox{ in } (0,NT)\times(0,1)\ .
\leqno(2.5) $$
Here $\,T\,$ is a positive constant.

\proclaim{Lemma 2.2}.
We put
$$
a(x)={x(1-x)\over 2}\,,
\ b(x)={1-2x\over 2}-NM(x)\,,
\ c(x)=-NM'(x)\,,
\leqno(2.6)$$
and define a differential operator $\,B\,$ by
$$
B={1\over2}a(x){\pd^2\over\pd x^2}+b(x){\pd\over\pd x}+c(x)\ .
\leqno(2.7)$$
If $\,v(t,x)\,$ is a solution of the equation
$$
{\pd v\over\pd t}=B v\hbox{ in } (0,T)\times(0,1)
\leqno(2.8)$$
then
$$
u(t,x)=\exp\bigl(-{\ds t\over2N}\bigr)\,v\bigl({\ds t\over N},x\bigr)
\leqno(2.9)$$
satisfies (2.5).

\Proof
Direct computation gives
$$
\eqalign{
Au
&={1\over2}{x(1-x)\over2N}\,u_{xx}
	+\Bigl({1-2x\over2N}-M(x)\Bigr)\,u_x
	-\Bigl({1\over2N}+M'(x)\Bigr)\,u\cr
&={1\over N}\Bigl({1\over2}a(x)\,u_{xx}+b(x)\,u_x+c(x)\,u-{1\over2}\,u\Bigr)\cr
&={1\over N}(Bu-{1\over2}\,u)\ .\cr
}
$$
Hence,
if $\,v(t,x)\,$ is a solution of (2.8)
, and if we put
$$
u(t,x)=\exp\bigl(-{\ds t\over2N}\bigr)\,v\bigl({\ds t\over N},x\bigr)\ ,
$$
then we have
$$
\eqalign{
u_t
&={1\over N}\exp\bigl(-{t\over2N}\bigr)\,v_t
	-{1\over2N}\exp\bigl(-{t\over2N})\,v\cr
&={1\over N}\exp\bigl(-{t\over2N}\bigr)\,Bv
	-{1\over2N}\exp\bigl(-{t\over2N})\,v\cr
&={1\over N}(Bu-{1\over2}\,u)\cr
&=Au\ .\qquad \QED\cr
}
$$

By virtue of Lemma 2.2,
the partial differential equation (2.5) is reduced to (2.8).
From view points of numerical analysis and population genetics,
it is much easier to solve (2.8) than (2.5).
In fact, the rate of fixation term, $\,\exp(-t/2N)\,$, is separated in (2.9)
and the original time scale is changed to a moderate one
in the process of reduction.

As will be shown later,
we can construct the desired {\it approximate general solution\/}
of the equation
$\,u_t=B u\,$
by using the following lemma recursively.

\proclaim{Lemma 2.3}.
Assume that
$\,u(t,x)\in C^{2,4}([0,\infty)\times[0,1])\,$
is a classical solution of the degenerate parabolic equation
$$
{\pd u\over\pd t}=B u\qquad (t>0,\ 0<x<1).
\leqno(2.10)$$
Then we obtain
$$
u(t,x)
=Lu(t,x)-h^4 Ru(t,x)
\leqno(2.11)$$
where
$\,0<h\ll 1\,$,
$\,0\leq x\pm\sqrt{a(x)}h\leq 1\,$,
$\,0\leq x+b(x)h^2\leq 1\,$,
and
$\,t-h^2\geq 0\,$.
Here
$$
\eqalign{
&Lu(t,x)\cr
&={1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)
	+{h^2\over3}\,c(x)u(t,x)\cr
}
\leqno(2.12)$$
and
$$
\eqalign{
&Ru(t,x)\cr
&={1\over3}\,
	\int_0^1\Bigl\{{a^2(x)\,(1-\theta )^3\over12}
	\Bigl(\,{\pd^4 u\over\pd x^4}\bigl(t,x+\sqrt{a(x)}\,\theta  h\bigr)
	+{\pd^4 u\over\pd x^4}\bigl(t,x-\sqrt{a(x)}\,\theta  h\bigr)\Bigr)\cr
&+b^2(x)\,(1-\theta )\,{\pd^2 u\over\pd x^2}(t,x+b(x)\,\theta  h^2)
	+(1-\theta )\,{\pd^2 u\over\pd t^2}(t-\theta  h^2,x)\Bigr\}\,d\theta \ .\cr
}
\leqno(2.13)$$

\Proof
(2.1) gives
$$
\eqalign{
&u\bigl(t,x\pm\sqrt{a(x)}h\bigr)\cr
&=u(t,x)\pm h\,a^{1/2}(x){\pd u\over\pd x}(t,x)\cr
&+{h^2\over2}\,a(x){\pd^2 u\over\pd x^2}(t,x)
	\pm{h^3\over6}\,a^{3/2}(x){\pd^3 u\over\pd x^3}(t,x)\cr
&+{h^4\over6}\,a^2(x)\int_0^1 (1-\theta )^3\,
	{\pd^4 u\over\pd x^4}\bigl(t,x\pm\theta \sqrt{a(x)}h\bigr)\,d\theta \ ;\cr
}
\leqno(2.14)$$
this implies
$$
\eqalign{
&{1\over2}u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over2}u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&=u(t,x)+{h^2\over2}\,a(x){\pd^2 u\over\pd x^2}(t,x)\cr
&+{h^4\over12}\,a^2(x)\int_0^1 (1-\theta )^3\,
	\Bigl(\,{\pd^4 u\over\pd x^4}\bigl(t,x+\theta \sqrt{a(x)}\,h\bigr)
	+{\pd^4 u\over\pd x^4}\bigl(t,x-\theta \sqrt{a(x)}\,h\bigr)\Bigr)
	\,d\theta \ .\cr
}
\leqno(2.15)$$
(2.1) also gives
$$
\eqalign{
&u(t,x+b(x)h^2)\cr
&=u(t,x)+h^2 b(x){\pd u\over\pd x}(t,x)\cr
&+h^4 b^2(x)\int_0^1 (1-\theta )\,
	{\pd^2 u\over\pd x^2}(t,x+\theta \,b(x)h^2)\,d\theta \cr
}
\leqno(2.16)$$
and
$$
\eqalign{
&u(t-h^2,x)\cr
&=u(t,x)-h^2{\pd u\over\pd t}(t,x)\cr
&+h^4\int_0^1 (1-\theta )\,
	{\pd^2 u\over\pd t^2}(t-\theta  h^2,x)\,d\theta \ .\cr
}
\leqno(2.17)$$
On the other hand,
$$
{1\over2}a(x){\pd^2 u\over\pd x^2}+b(x){\pd u\over\pd x}
	-{\pd u\over\pd t}=-c(x)u
\leqno(2.18)$$
follows from (2.10).
Combining (2.18) with (2.15), (2.16) and (2.17), we obtain
$$
\eqalign{
&{1\over2}u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
+{1\over2}u\bigl(t,x-\sqrt{a(x)}\,h\bigr)
	+u(t,x+b(x)h^2)
	+u(t-h^2,x)\cr
&=3u(t,x)-h^2c(x)u(t,x)\cr
&+h^4\Biggl\{
	{a^2(x)\over12}\int_0^1 (1-\theta )^3\,
	\Bigl(\,{\pd^4 u\over\pd x^4}\bigl(t,x+\theta \sqrt{a(x)}\,h\bigr)
	+{\pd^4 u\over\pd x^4}\bigl(t,x-\theta \sqrt{a(x)}\,h\bigr)\Bigr)
	\,d\theta \cr
&+b^2(x)\int_0^1 (1-\theta )\,
	{\pd^2 u\over\pd x^2}(t,x+\theta \,b(x)h^2)\,d\theta 
	+\int_0^1 (1-\theta )\,
	{\pd^2 u\over\pd t^2}(t-\theta  h^2,x)\,d\theta \Biggr\}\ .\cr
}
\leqno(2.19)$$
This completes the proof.
\QED

Lemma 2.3\ \ shows that if $\,u(t,x)\,$ is a $C^{2,4}$ smooth solution,
then $\,u(t,x)\,$ can be rewritten as
$$
\eqalign{
u(t,x)
&=Lu(t,x)+O(h^4)\cr
&={1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)
	+{h^2\over3}\,c(x)u(t,x)
	+O(h^4)\ ;\cr
}
\leqno(2.20)$$
recursive use of this formula gives the desired {\it approximate general
solution\/} of degenerate parabolic equation (2.10).
Here, by (2.10) and (2.12), $\,O(h^4)\,$ depends on
$\,u_{xxxx},\ u_{xx}\,$ and $\,B^2u\,$
and further,
a standard method of maximum principle ensures that we can estimate
$\ \pd_x^ku\ \ (0\leq  k\leq  4)\ $
explicitly (cf.\/, for example, [5]).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%					%
%		Section 3		%
%					%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigbreak
\centerline{\bf \S3. Construction of Approximate General Solution}
\medskip\nobreak

In this section,
we construct approximate general solution of (2.10)
and give error bounds (Theorems 3.1 and 3.2).
We also show that our result can be translated
into symbolic list operations quite easily (Reduce Program).

We note that because of the phenomenon which (1.1) describes,
no boundary condition should be given on $\ x=0,\ 1\ $,
i.\ e., the boundary $\ x=0,\ 1\ $ should be either the {\it entrance\/}
or {\it natural boundary\/} of a diffusion process generated by (2.8).
So, without loss of practical generality, we may assume that
$$
b(x)\cases{
\geq 0 & in a neighborhood of 0 in [0,1],\cr
\cr
\leq 0 &  in a neighborhood of 1 in [0,1].\cr
}
\leqno(3.1)$$
Furthermore, we may consider
$$
c(x)\leq 0\quad in\ [0,1]
\leqno(3.2)$$
since if we put
$$
w(t,x)=e^{-tc_0}u(t,x)\ ,\quad c_0=\sup_{0\leq  x\leq 1}c(x)
$$
then
$$
u_t=Bu\ \Longleftrightarrow\ w_t=(B-c_0)w\ .
$$
Throughout this section,
$\,u(t,x)\in C^{2,4}([0,\infty)\times[0,1])\,$
is a classical solution of the equation (2.10) and
$\,h\,$ denotes a sufficiently small fixed positive constant,
i.\ e., $\,0<h\ll1\,$.

By using (2.20), if we replace
$\,(h^2/3)c(x)u(t,x)\,$
with
$$
\eqalign{
&{h^2\over3}c(x)\Bigl(\,{1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)
	+{h^2\over3}\,c(x)u(t,x)
	+O(h^4)\Bigr)\cr
}
$$
in (2.20), we obtain
$$
\eqalign{
u(t,x)
&=\Bigl(1+{h^2\over3}c(x)\Bigr)
	\Bigl(\,{1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)\,\Bigr)
	+{h^4\over9}\,c^2(x)u(t,x)
	+O(h^4)\ .\cr
}
\leqno(3.3)$$
Here, since
$\,c^2(x)|u(t,x)|=\bigl(NM'(x)\bigr)^2|u(t,x)|\,$
may not be so small enough in some cases, it is better not to replace
$\,(h^4/9)\,c^2(x)u(t,x)\,$
with $\,O(h^4)\,$.
So, we repeat the above recursive transformation again for (3.3),
i.\ e., we substitute
$$
\eqalign{
{h^4\over9}\,c^2(x)\Bigl\{
&\Bigl(1+{h^2\over3}c(x)\Bigr)
	\Bigl(\,{1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)\,\Bigr)
	+{h^4\over9}\,c^2(x)u(t,x)
	+O(h^4)\Bigr\}\cr
}
$$
for
$\ \ds{h^4\over9}\,c^2(x)u(t,x)\ $
in the equation (3.3),
and get the following:
$$
\eqalign{
u(t,x)
&=\Bigl(1+{h^2\over3}c(x)+{h^4\over9}c^2(x)\Bigr)
	\Bigl(\,{1\over6}\,u\bigl(t,x+\sqrt{a(x)}\,h\bigr)
	+{1\over6}\,u\bigl(t,x-\sqrt{a(x)}\,h\bigr)\cr
&+{1\over3}\,u\bigl(t,x+b(x)h^2\bigr)
	+{1\over3}\,u(t-h^2,x)\,\Bigr)
	+O(h^4)\ .\cr
}
\leqno(3.4)$$
This is the key formula of this section.

We define a function $\ 0<\varepsilon(x)\leq  h\ $ by
$$
\varepsilon(x)=\cases{
\quad h &if\quad $h^2\leq  x\leq  1-h^2$
	\enspace or\enspace $x=0,\ 1$\cr
\cr
\ds{x\over\sqrt{a(x)}} &if\quad $0<x<h^2$\cr
\cr
\cr
\ds{1-x\over\sqrt{a(x)}} &if\quad $1-h^2<x<1$\ .\cr
}
\leqno(3.5)$$
Here it is to be noted that, by (2.6), (3.1) and (3.5),
$$
x\pm\sqrt{a(x)}\,\varepsilon(x)\in [0,1]\hbox{ and } x+b(x)\,\varepsilon^2(x)\in [0,1]
\leqno(3.6)$$
holds for any $\,0\leq  x\leq  1\,$.

\proclaim{Lemma 3.1}.
For $\,0<h\ll 1\,$,
$$
h^2\leq  x\leq  1-h^2	\hbox{ implies }
	{h^2\over4}\leq  x\pm\sqrt{a(x)}\,h\leq  1-{h^2\over4}\ .
\leqno(3.7)$$

\Proof
We put
$$
f(x)=x-\sqrt{a(x)}\,h-{h^2\over 4}=x-\sqrt{x(1-x)\over 2}\,h-{h^2\over 4}\ .
$$
Direct computation gives
$$
f'(x)>0\hbox{ for } {1\over 2}-{1\over\sqrt{2h^2+4}}<x\leq  1
$$
and
$$
f(h^2)>0\ ,\quad h^2>{1\over 2}-{1\over\sqrt{2h^2+4}}\ .
$$
Hence, we obtain
$$
f(x)>0\hbox{ for } h^2\leq  x\leq  1\ ;
$$
this implies
$$
x\pm\sqrt{a(x)}\,h\geq {h^2\over 4}\hbox{ for } h^2\leq  x\leq  1\ .
$$
Similarly, we obtain
$$
x\pm\sqrt{a(x)}\,h\leq 1-{h^2\over 4}\hbox{ for } 0\leq  x\leq  1-h^2\ .\quad\QED
$$

\proclaim{Lemma 3.2}.
For $\,0<h\ll 1\,$,
$$
{h^2\over4}\leq  y\leq  1-{h^2\over4}
	\hbox{ implies }
	\varepsilon(y)\geq {h\over\sqrt{2}}
\leqno(3.8)$$

\Proof
Seeing the definition of $\varepsilon(x)$,
we have only to show (3.8) for
$\,\ds{h^2\over4}\leq  y<h^2\,$ and $\,\ds 1-h^2<y\leq  1-{h^2\over4}\ $.
We put
$$
g(y)=\varepsilon^2(y)={2y\over1-y}\hbox{ for }{h^2\over 4}\leq y\leq h^2\ .
$$
Since $\,g\kern0.07em'(y)>0\,$ for $\,\ds{h^2\over4}\leq  h\leq  h^2\,$,
we have
$$
\min_{ h^2/4\leq  y\leq  h^2}g(y)=g\Bigl({h^2\over4}\Bigr)
	={2h^2\over4-h^2}\geq {h^2\over2}\ ;
$$
this gives
$$
\varepsilon(y)\geq {h\over\sqrt{2}}\hbox{ for }{h^2\over4}\leq  y<h^2\ .
$$
Similarly, we obtain
$$
\varepsilon(y)\geq {h\over\sqrt{2}}\hbox{ for }1-h^2<y\leq {1-{h^2\over4}}\ .\quad\QED
$$

By $\,M_{\varepsilon}\,$, we denote a difference operator
$$
M_{\varepsilon}f(t,x)=\cases{
\ds\Bigl(1+{\varepsilon^2(x)\over3}c(x)+{\varepsilon^4(x)\over9}c^2(x)\Bigr)\cr
\ds\times\Bigl(\,{1\over6}\,f\bigl(t,x+\sqrt{a(x)}\,\varepsilon(x)\bigr)
	+{1\over6}\,f\bigl(t,x-\sqrt{a(x)}\,\varepsilon(x)\bigr)\cr
\ds\ +{1\over3}\,f\bigl(t,x+b(x)\varepsilon^2(x)\bigr)
	+{1\over3}\,f(t-\varepsilon^2(x),x)\,\Bigr)
	&if\enspace $t\geq h^2$\cr
\cr
\cr
\ds\ f(t,x)
	&otherwise\ .\cr
}
\leqno(3.9)$$
We define functions $\,p_k(t,x;s,y),\ k=0,1,2,\cdots\,$ as follows :
$$
p_0(t,x,s,y)=\cases{
\ 1 &if\enspace $(s,y)=(t,x)$\cr
\cr
\ 0 &otherwise\ ,\cr
}
\leqno(3.10)$$
$$
p_1(t,x,s,y)=\cases{
\ds{1\over6}\Bigl(1+{\varepsilon^2(x)\over3}c(x)+{\varepsilon^4(x)\over9}c^2(x)\Bigr)
	&if\enspace $(s.y)=\big(t,x\pm\sqrt{a(x)}\,\varepsilon(x)\big)$\cr
\cr
\ds{1\over3}\Bigl(1+{\varepsilon^2(x)\over3}c(x)+{\varepsilon^4(x)\over9}c^2(x)\Bigr)
	&if\enspace $(s.y)=\big(t,x+b(x)\,\varepsilon^2(x)\big)$\cr
\cr
\ds{1\over3}\Bigl(1+{\varepsilon^2(x)\over3}c(x)+{\varepsilon^4(x)\over9}c^2(x)\Bigr)
	&if\enspace $(s.y)=\big(t-\varepsilon^2(x),x\big)$\cr
\cr
\ 0 &otherwise\ ,\cr
}
\leqno(3.11)$$
and $\ p_{k+1}(t,x,s,y)\ \ (k\geq 1)\ $ is a function such that
$$
\int_D f(s,y)\,p_{k+1}(t,x,ds,dy)
	=\int_D M_{\varepsilon}f(s,y)\,p_k(t,x,ds,dy)
\leqno(3.12)$$
$$
\Biggl(\quad i.\ e.,\quad\sum_{(s,y)\in D}f(s,y)\,p_{k+1}(t,x,s,y)
	=\sum_{(s,y)\in D}M_{\varepsilon}f(s,y)\,p_k(t,x,s,y)\quad\Biggr)
$$
for any function $\,f(t,x)\,$ defined in $D\equiv [0,\infty]\times[0,1]\,$.
Here (3.12) is well-defined, since
$\,\sqrt{a(x)}\,\varepsilon(x)\not=|b(x)|\,\varepsilon^2(x)\,$
for $\,0<h\ll1$. It is clear that $\,p_{k+1}(t,x,s,y)\, $ in (3.12) is 
uniquely determined by $\,p_k(t,x,s,y)\,$.

Then (3.10) implies
$$
u(t,x)=\int_D p_0(t,x,ds,dy)\,u(s,y)\ .
\leqno(3.13)$$
Since $\,\varepsilon^4(x)=O(h^4)\,$, (3.4) and (3.11) give
$$
u(t,x)=\int_D p_1(t,x,ds,dy)\,u(s,y)+O(h^4)\ .
\leqno(3.14)$$
Furthermore, by (3.4), (3.12) and (3.14), we have
$$
\eqalign{
u(t,x)
&=\int_D p_1(t,x,ds,dy)\,u(s,y)+O(h^4)\cr
&=\int_D p_1(t,x,ds,dy)\big(M_{\varepsilon}u(s,y)+O(h^4)\big)+O(h^4)\cr
&=\int_D p_2(t,x,ds,dy)\,u(s,y)+O(h^4)r+O(h^4)\cr
&=\int_D p_2(t,x,ds,dy)\big(M_{\varepsilon}u(s,y)+O(h^4)\big)+O(h^4)r+O(h^4)\cr
&=\int_D p_3(t,x,ds,dy)\,u(s,y)+O(h^4)r^2+O(h^4)r+O(h^4)\ ,\cr
}
\leqno(3.15)$$
where
$$
r=\sup_{0\leq  x\leq 1}
	\Bigl(1+{\varepsilon^2(x)\over3}c(x)+{\varepsilon^4(x)\over9}c^2(x)\Bigr)
	=1+O(h^4)
\leqno(3.16)$$
by virtue of (3.2) and (3.5).
We can repeat this procedure inductively.

Combining the above results, we obtain the following

\proclaim{Theorem 3.1}.
Assume that $\,u(t,x)\in C^{2,4}([0,\infty)\times[0,1])\,$
is a classical solution of the degenerate parabolic equation (2.10)
and $\,0<h\ll 1\,$ is a fixed small positive constant.
Then,
for $\,(t,x)\in (0,\infty)\times[0,1]\,$,
we have
$$
u(t,x)=\int_{0\leq s\leq t,\ 0\leq y\leq 1}p_k(t,x,ds,dy)\,u(s,y)
	+O(h^4)\sum_{\nu =0}^{k-1}r^{\nu }
\leqno(3.17)$$
for any $\,k=0,1,2,\cdots$\ .
Here $O(h^4)$ is independent of $\,k\,$.

In (3.17), if we replace $\,u(s,y),\ 0\leq  s<h^2\,$
with given initial data $\,\phi(y)\,$,
we get the following

\proclaim{Theorem 3.2}.
Assume that $\,u(t,x)\in C^{2,4}([0,\infty)\times[0,1])\,$
is a classical solution
of  the degenerate parabolic equation (2.10) satisfying initial condition
$$
u(0,x)=\phi(x)\qquad (0<x<1)\,,
$$
and $\,0<h\ll 1\,$ is a fixed small positive constant.
Then,
for $\,(t,x)\in (0,\infty)\times[0,1]\,$,
we have
$$
\eqalign{
u(t,x)
&=\int_{0\leq  s<h^2,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\,\phi(y)\cr
&+O(h^2)\,r^k+O(1)\,r^k\sum_{\ell=0}^{k_0}{\,k\,\choose\,\ell\,}
	\Bigl({1\over3}\Bigr)^\ell\Bigl({2\over3}\Bigr)^{k-\ell}
	+O(h^4)\sum_{\nu =0}^{k-1}r^{\nu }\cr
}
\leqno(3.18)$$
for any $\,k\geq 2t/h^2\,$,
where $k_0$ is the smallest integer satisfying
$\,2t/h^2\leq  k_0\leq  k\,$.
Here $O(1),\ O(h^2)$ and $O(h^4)$ are independent of $\,k\,$.

Hence, we obtain the desired {\it approximate general solution\/}
$$
u(t,x)\sim\int_{0\leq  s<h^2,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\,\phi(y)
\leqno(3.19)$$
for $\,k\gg 2t/h^2\,$.
For a given $\,0<h\ll 1\,$ and $\,(t,x)\in (0,\infty)\times[0,1]\,$,
it will suffice to choose a certain $\,k\gg 2t/h^2\,$ which minimizes
the error
$$
h^2r^k+r^k\sum_{\ell=0}^{k_0}{\,k\,\choose\,\ell\,}
	\Bigl({1\over3}\Bigr)^\ell\Bigl({2\over3}\Bigr)^{k-\ell}
	+h^4\sum_{\nu =0}^{k-1}r^{\nu }\ .
$$
So, computing the above error for $\,2t/h^2\leq  k\leq  1/h^4\,$,
we are able to know when we have to stop the iteration in advance.

\smallbreak\noindent{\it Proof of Theorem 3.2\/.\enspace}
Theorem 3.1 gives
$$
\eqalign{
u(t,x)
&=\int_{0\leq  s<h^2,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\,\phi(y)+O(h^2)\,r^k\cr
&+\int_{h^2\leq  s\leq 1,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\,u(s,y)
	+O(h^4)\sum_{\nu =0}^{k-1}r^{\nu }\ .\cr
}
\leqno(3.20)$$
Seeing the definition of $\,p_k(t,x,s,y)\,$,
we obtain, by Lemmas 3.1 and 3.2,
$$
\eqalign{
&\Bigl|\,\int_{h^2\leq  s\leq 1,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\,u(s,y)\,\Bigr|\cr
&\leq  C_1\int_{h^2\leq  s\leq 1,\ 0\leq  y\leq 1}p_k(t,x,ds,dy)\cr
&\leq  C_2\,r^k\sum_{\ell=0}^{k_0}{\,k\,\choose\,\ell\,}
	\Bigl({1\over3}\Bigr)^\ell\Bigl({2\over3}\Bigr)^{k-\ell}\cr
}
$$
where $\ C_i\ (i=1,2)\ $ are positive constants independent of $\,k\,$.
\QED

If we identify a linear combination of the form
$$
\int_{0\leq s\leq t,\ 0\leq y\leq 1}p_k(t,x,ds,dy)\,u(s,y)=
	q_1 u(t_1, x_1)+q_2 u(t_2, x_2)+\cdots+ q_n u(t_n, x_n)
$$
with a list
$$
\bigl((q_1\ t_1\ x_1)\ (q_2\ t_2\ x_2)\ \cdots\ (q_n\ t_n\ x_n)\bigr)
$$
of LISP type,
then the essential part of Theorem 3.2 is translated into the following
list operations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% beginning of REDUCE program %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent
{\bf Reduce Program.}

\ifx\gtfam\undefined
  \ifx\dm\undefined
    \ifx\tendm\undefined
      \def\mc{\null}
    \else
      \def\mc{\tendm}
    \fi
  \else
    \def\mc{\dm}
  \fi
  \ifx\dg\undefined
    \ifx\tendg\undefined
      \def\gt{\null}
    \else
      \def\gt{\tendg}
    \fi
  \else
    \def\gt{\dg}
  \fi
\fi
\ifx\sc\undefined
  \def\sc{\null}
\fi

\tt\mc 

\noindent
p\kern-.45em p\kern-.05em r\kern-.45em r\kern-.05em o\kern-.45em 
o\kern-.05em c\kern-.45em c\kern-.05em e\kern-.435em e\kern-.065em 
d\kern-.45em d\kern-.05em u\kern-.45em u\kern-.05em r\kern-.45em 
r\kern-.05em e\kern-.435em e\kern-.065em {\tt\mc \kern0.500em}hybrid{
\tt\_\kern.141em}method(t,{\tt\mc \kern0.500em}x,{\tt\mc \kern0.500em}n);

\noindent
b\kern-.45em b\kern-.05em e\kern-.435em e\kern-.065em g\kern-.45em 
g\kern-.05em i\kern-.45em i\kern-.05em n\kern-.45em n\kern-.05em ;

\noindent
{\tt\mc \kern1.500em}list{\tt\_\kern.141em}in{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}\rm\mc $\kern.024em\{$$\kern.024em\{$1,{\tt\mc 
\kern0.500em}t,{\tt\mc \kern0.500em}x$\}\kern.024em$\tt\mc $\}\kern.024em$;

\noindent
{\tt\mc \kern1.500em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}\rm\mc $\kern.024em\{$$\}\kern.024em$\tt\mc ;

\noindent
{\tt\mc \kern1.500em}w\kern-.45em w\kern-.05em h\kern-.45em h\kern-.05em 
i\kern-.45em i\kern-.05em l\kern-.45em l\kern-.05em e\kern-.435em 
e\kern-.065em {\tt\mc \kern0.500em}n{\tt\mc \kern0.500em}{\tt >}{\tt\mc 
\kern0.500em}0{\tt\mc \kern0.500em}d\kern-.45em d\kern-.05em o\kern-.45em o\kern-.05em 

\noindent
{\tt\mc \kern1.500em}{\tt <}{\tt <}

\noindent
{\tt\mc \kern3.000em}w\kern-.45em w\kern-.05em h\kern-.45em h\kern-.05em 
i\kern-.45em i\kern-.05em l\kern-.45em l\kern-.05em e\kern-.435em 
e\kern-.065em {\tt\mc \kern0.500em}length(list{\tt\_\kern.141em}in)
{\tt\mc \kern0.500em}{\tt >}{\tt\mc \kern0.500em}0{\tt\mc 
\kern0.500em}d\kern-.45em d\kern-.05em o\kern-.45em o\kern-.05em 

\noindent
{\tt\mc \kern3.000em}{\tt <}{\tt <}

\noindent
{\tt\mc \kern4.500em}tmp{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}first(list{\tt\_\kern.141em}in);

\noindent
{\tt\mc \kern4.500em}p{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}first(tmp);

\noindent
{\tt\mc \kern4.500em}s{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}first(rest(tmp));

\noindent
{\tt\mc \kern4.500em}y{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}first(rest(rest(tmp)));

\noindent
{\tt\mc \kern4.500em}q{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}(1+h{\tt *}{\tt *}2{\tt *}c(y){\tt /}3+h{\tt *}{\tt 
*}4{\tt *}c(y){\tt *}{\tt *}2{\tt /}9){\tt *}p;

\noindent
{\tt\mc \kern4.500em}i\kern-.45em i\kern-.05em f\kern-.45em 
f\kern-.05em {\tt\mc \kern0.500em}domain{\tt\_\kern.141em}p(s){\tt\mc 
\kern0.500em}neq{\tt\mc \kern0.500em}0{\tt\mc \kern0.500em}t\kern-.44em 
t\kern-.06em h\kern-.45em h\kern-.05em e\kern-.435em e\kern-.065em 
n\kern-.45em n\kern-.05em 

\noindent
{\tt\mc \kern4.500em}{\tt <}{\tt <}

\noindent
{\tt\mc \kern6.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}cons(\rm\mc $\kern.024em\{$q{\tt /}6,{\tt\mc 
\kern0.500em}s,{\tt\mc \kern0.500em}y+sqrt(a(y)){\tt 
*}h$\}\kern.024em$\tt\mc ,{\tt\mc \kern0.500em}list{\tt\_\kern.141em}tmp);

\noindent
{\tt\mc \kern6.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}cons(\rm\mc $\kern.024em\{$q{\tt /}6,{\tt\mc 
\kern0.500em}s,{\tt\mc \kern0.500em}y{\tt -}sqrt(a(y)){\tt *}h$\}
\kern.024em$\tt\mc ,{\tt\mc \kern0.500em}list{\tt\_\kern.141em}tmp);

\noindent
{\tt\mc \kern6.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}cons(\rm\mc $\kern.024em\{$q{\tt /}3,{\tt\mc 
\kern0.500em}s,{\tt\mc \kern0.500em}y+b(y){\tt *}h{\tt *}{\tt *}2$\}
\kern.024em$\tt\mc ,{\tt\mc \kern0.500em}list{\tt\_\kern.141em}tmp);

\noindent
{\tt\mc \kern6.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}cons(\rm\mc $\kern.024em\{$q{\tt /}3,{\tt\mc 
\kern0.500em}s{\tt -}h{\tt *}{\tt *}2,{\tt\mc \kern0.500em}y$\}
\kern.024em$\tt\mc ,{\tt\mc \kern0.500em}list{\tt\_\kern.141em}tmp)

\noindent
{\tt\mc \kern4.500em}{\tt >}{\tt >}

\noindent
{\tt\mc \kern4.500em}e\kern-.435em e\kern-.065em l\kern-.45em l\kern-.05em
 s\kern-.45em s\kern-.05em e\kern-.435em e\kern-.065em 

\noindent
{\tt\mc \kern6.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}cons(\rm\mc $\kern.024em\{$p,{\tt\mc \kern0.500em}s,
{\tt\mc \kern0.500em}y$\}\kern.024em$\tt\mc ,{\tt\mc 
\kern0.500em}list{\tt\_\kern.141em}tmp);

\noindent
{\tt\mc \kern4.500em}list{\tt\_\kern.141em}in{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}rest(list{\tt\_\kern.141em}in)

\noindent
{\tt\mc \kern3.000em}{\tt >}{\tt >};

\noindent
{\tt\mc \kern3.000em}list{\tt\_\kern.141em}in{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}list{\tt\_\kern.141em}tmp;

\noindent
{\tt\mc \kern3.000em}list{\tt\_\kern.141em}tmp{\tt\mc \kern0.500em}:=
{\tt\mc \kern0.500em}\rm\mc $\kern.024em\{$$\}\kern.024em$\tt\mc ;

\noindent
{\tt\mc \kern3.000em}n{\tt\mc \kern0.500em}:={\tt\mc 
\kern0.500em}n{\tt -}1

\noindent
{\tt\mc \kern1.500em}{\tt >}{\tt >};

\noindent
{\tt\mc \kern1.500em}return{\tt\mc \kern0.500em}list{\tt\_\kern.141em}in

\noindent
e\kern-.435em e\kern-.065em n\kern-.45em n\kern-.05em d\kern-.45em 
d\kern-.05em ;

\noindent
\rm\mc

%%%%%%%%%%%%%%%%%%%%%%%%%
% end of REDUCE program %
%%%%%%%%%%%%%%%%%%%%%%%%%

Here {\tt\ domain\_p(t)\ } is a function
defined in $\,[0, \infty)\,$ such that
$$
{\tt domain\_p(t)}=\cases{
	\ 1 &if\enspace$t\geq h^2$\cr
	\cr
	\ 0 &otherwise\ .\cr}
$$
Functions $\,a(x),\ b(x)\,$ and $\,c(x)\,$ are coefficients defined by 
(2.6).
To be more explicit, the parameter $\,h\,$ should be modified
in each step of main loop of the above Reduce Program so that
$$
y\pm\sqrt{a(y)}\,h,\ y+b(y)\,h^2\in [0,1]
$$
is valid and also, we should identify any pair of list
$\,(q_i\ t_i\ x_i)\,$ and $\,(q_j\ t_j\ x_j)\,$
with
$\,(q_i+q_j\ t_i\ x_i)\,$
when $\,(t_i\ x_i)=(t_j\ x_j)\,$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%					%
%		Section 4		%
%					%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigbreak
\centerline{\bf \S4. Examples}\medskip\nobreak

We shall solve the following two problems by using (3.19) :
$$
\cases{
\ds{\pd u\over\pd t}={1\over4N}{\pd^2\over\pd x^2}\big(x(1-x)u\big)
	\qquad (0<t\leq 4N,\ 0<x<1)\cr
\cr
u(0,x)\sim\delta(x-0.5)\ ,\cr
}
\leqno(4.1)$$
$$
\cases{
\ds{\pd u\over\pd t}={1\over4N}{\pd^2\over\pd x^2}\big(x(1-x)u\big)
	\qquad (0<t\leq 4N,\ 0<x<1)\cr
\cr
u(0,x)\sim\delta(x-0.1)\ .\cr
}
\leqno(4.2)$$
These problems describe the simplest situation
in which a pair of alleles $\alpha _1$ and $\alpha _2$ are segregating
with frequences $\,x\,$ and $\,1-x\,$ respectively
in a randomly mating population of N monoecious individuals,
and in which only the random sampling of gametes in reproduction
causes gene frequency change.

Let $\,u(t,x)\,$ be a classical solution of (4.1)
and let $\,u_k(t,x)=\sum_j q_j\phi(x_j)\,$ be
an approximate solution of (4.1) constructed
for sufficiently small $\,h>0\,$.
Then, since $\,c(x)=-NM'(x)=0\,$ implies $\,r=1\,$,
Theorem 3.2 gives
$$
|u(t,x)-u_k(t.x)|\leq  
C\,\Biggl(h^2+\sum_{\ell=0}^{[2t/h^2]+1}{\,k\,\choose\,\ell\,}
	\Bigl({1\over3}\Bigr)^\ell\Bigl({2\over3}\Bigr)^{k-\ell}
	+h^4k\Biggr)
\leqno(4.3)$$
for $\,k\geq 2t/h^2\,$.
Here we note that, seeing the proof of Lemma 2.3,
we can estimate the above constant $\,C\,$ by using maximum principle
and also, according to the fundamental property of binomial distribution,
the second term of the right hand side of (4.3) tends to 0
when $\,k\to\infty\,$.
Thus, if we choose $\,2t/h^2\ll k\ll 1/h^4\,$ appropriately,
$\,|u(t,x)-u_k(t.x)|\,$ would become sufficiently small.
The same argument remains true for the problem (4.2).

The graphs of our approximate solutions of (4.1) and (4.2) are given
in Figures 1 and 2.
It is to be noted that our solutions are approximately equal to
Kimura's exact ones
and they lead to the same conclusion that Kimura obtained
\ (cf.\ [$\,$2, Chapter 8$\,$]).

%%%%%%%%%%%%
% Figure 1 %
%%%%%%%%%%%%

\topinsert\vskip 8cm
\special{psfile=sol1.eps hscale=70 vscale=70}
\centerline{Figure 1} \endinsert


%%%%%%%%%%%%
% Figure 2 %
%%%%%%%%%%%%

\topinsert\vskip 11.5cm
\special{psfile=sol2.eps hscale=60 vscale=60}
\centerline{Figure 2} \endinsert


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%						%
%		Acknowledgment			%
%						%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\noindent{\bf Acknowledgement}.
The author would like to express his hearty thanks to the referee
for his helpful criticisms, suggestions and encouragement.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%					%
%		References		%
%					%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vfill\eject
\centerline{\bf References}\medskip\nobreak
\item{[1]}
K. Amano,
Numerical-symbolic hybrid method for biharmonic Dirichlet problem,
({\it to appear\/}).

\item{[2]}
J. F. Crow and M. Kimura,
An Introduction to Population Genetics Theory,
{\sl Burgess Publishing Company}, 1970.

\item{[3]}
W. J. Ewens,
Mathematical Population Genetics,
{\sl Springer-Verlag}, 1989.

\item{[4]}
A. C. Hearn,
REDUCE User's Manual, version 3.5,
{\sl RAND Publication}, CP78 (Rev. 7/94), 1994.

\item{[5]}
O. A. Oleinik and E. V. Radkevich,
Second Order Equations with Nonnegative Characteristic Form,
{\sl Amer. Math. Soc., Province, Rhode Island and Plenum Press,
New York}, 1973.



{\sc Department of Mathematics, Josai University, Sakado, Saitama, JAPAN}

E-mail: kamano@tansei.cc.u-tokyo.ac.jp

\bye