
%%%MUST have jsc.sty in the same directory, also obtainable in Zeilberger's
%%%public directory
%%%%LaTeX file for a paper by C.Orr and D.Zeilberger, to appear
%%in Journal of Symbolic Computation

\documentstyle{jsc}
\title[A Computer Algebra Approach to the Discrete Dirichlet Problem]{A 
Computer Algebra Approach to the Discrete Dirichlet Problem}
\author[Craig R. Orr and Doron Zeilberger]{CRAIG R. ORR AND DORON ZEILBERGER\\
Department of Mathematics, Temple University, Philadelphia, PA 19122, USA}
\begin{document}
\maketitle

\begin{abstract}
The Gr\"obner Basis Algorithm is used to find closed form expressions
for the generating functions of finite difference equations. Such difference
equations arise in elliptic PDE's, random walk problems, and gambler's 
ruin problems.
\end{abstract}


\section*{Introduction}
Solutions to the Dirichlet problem for (elliptic) partial difference
equations are needed
both in the solutions of boundary value problems of elliptic
partial differential equations (Birkhoff,1972) and in random walk problems
(Doyle and Snell, 1984).

In order to approximate solutions of boundary value problems in partial
differential equations,
one sets up a so-called finite difference scheme on a discrete grid
and considers the resulting finite difference problem. The solution
gives approximations for the values of the desired function at the
grid points. The finer the grid, the better the approximation. However,
the effort also grows, since the number of mesh points grows.

In random walk problems, the value of interest is either the probability
of getting from a fixed point to another point under a constrained walk,
or the expected life, as in gambler's ruin problems. The solution can
be formulated in terms of an appropriate boundary value problem for
a certain region of the lattice. 

As far as we know, the only method  that is used for obtaining exact solutions
for such
discrete Dirichlet problems is by straightforward linear algebra, getting
a large, but sparse, system of linear equations. A whole branch of numerical
linear algebra that treats sparse systems was developed from the need
to solve such boundary value problems.

Here we propose a new method for solving such problems, that reduces
the problem to a much smaller system of linear equations ( $O(N)$
as opposed to $O(N^{2})$), albeit no longer sparse.

Our method uses a very small part of Buchberger's celebrated method
of Gr\" obner basis, namely taking the ``normal form'' of a polynomial
(Buchberger,1976). This method
is implemented, among other places, in {\sc maple} in the {\tt grobner}
package.


We will illustrate our method by a detailed example:
The two-dimensional gambler's ruin problem, for which no closed form solution
is known to exist, and probably doesn't exist. 
The manipulations can be mechanized, but we put
in the details for expository sake.

\section*{Gambler's Ruin Examples \\ with Generating Functions}
Consider a gambler who wins a dollar or loses a dollar with probability 
$\frac{1}{2}$. The gambler starts out with $i$ dollars and an opponent starts
out with $n-i$ dollars. The game ends when one of the two players has no money 
left. How long should the game be expected to last? This problem is known as 
the symmetric classical gambler's ruin problem. The expected duration is known
and is given by $f_{n}(i)=i (n-i)$ (Feller, 1950, pp 311--318).

Now consider a two-dimensional version of this game: the gambler has two piles 
of money, say $i$ dollars
and $j$ shekels. The opponent has $M-i$ dollars and $N-j$ shekels. In this 
game, the gambler wins a dollar, loses a dollar, wins a shekel, or loses a 
shekel each with probabiliy $\frac{1}{4}$. The game ends when either the 
gambler or the opponent has one pile reduced to nothing. How long should this 
game be expected to last? 

Consider first the one-dimensional(pile) game. Let 
$f(i)= f_{n}(i)$ be the expected duration of the game. Then 
$$2 f(i) - f(i+1) - f(i-1)=2.$$
Let $$F(x)=\sum _{i=1}^{n-1}f(i)x^{i}$$ 
be the generating function for the problem. Then,

$$(2-x-x^{-1})F(x)  = 2\sum _{i=1}^{n-1}x^{i}-f(n-1)x^{n}-f(1).$$

Now let $f(n-1)=f(1)=c$. Then

$$(2-x-x^{-1})F(x)=\frac{2(x-x^{n})}{1-x} - c(x^{n}+1).$$
So,
$$F(x) ={\frac {\left (2\,x-2\,x^{n}-cx^{n}+cx^{n+1}-c+cx\right )x}{\left (x-1
\right )^{3}}}.$$

Now $F$, we know, is a polynomial by definition. So it can be determined
that $c=n-1$ by long division and setting the denominator equal to $0$.
( Alternatively, we could use l'Hopital's rule to find that
$c=n-1$. ) Therefore,
$$F(x) = {\frac {x\left (2\,x-2\,x^{n}-\left (n-1\right )x^{n}+\left (n-1
\right )x^{n+1}-n+1+\left (n-1\right )x\right )}{\left (x-1\right )^{3
}}}.$$
Therefore, we can find $F(i)$ by finding the constant term of 
$F(x)/x^{i}$.

\begin{eqnarray*}
f(i) & = & CT \left ( -{\frac {\left (n+1\right )x^{2}-\left (n-1\right )x}{\left (1-x
\right )^{3}x^{i}}} \right ) \\
& = & (n-i)i,
\end{eqnarray*}
which is the expected result.

The two-dimensional(pile) case is more difficult, since we will not be able to 
do simple long division at the end to get the result. 

Let $f(i,j)=f_{M,N}(i,j)$ be the expected duration of the game when the gambler
has $i$ dollars and $j$ shekels. Then
$$4f(i,j) - f(i+1,j)-f(i-1,j)-f(i,j+1)-f(i,j-1)=4$$ with 
$$f(i,0)=f(i,N)=f(0,j)=f(M,j)=0.$$
Let $$F(x,y) = F_{M,N}(x,y) = 
\sum_{i=1}^{M-1} \sum_{j=1}^{N-1} f(i,j) x^{i}y^{j}.$$ 
Then,

\begin{eqnarray*}
(4-x-x^{-1}-y-y^{-1}) F(x,y) 
&=&{\frac {\left (4\,x-4\,x^{M}\right )\left (y-y^{N}\right )}{\left (x-1
\right )\left (y-1\right )}} - x^{M}\sum_{j=1}^{N-1}f(M-1,j)y^{j} - \\
& &\sum_{j=1}^{N-1}f(1,j)y^{j} - 
y^{N}\sum_{i=1}^{M-1}f(i,N-1)x^{i} - \sum_{i=1}^{M-1}f(i,1)x^{i}.
\end{eqnarray*}
Now let $$\phi (x)= \sum_{i=1}^{M-1}f(i,N-1)x^{i}= \sum_{i=1}^{M-1}f(i,1)x^{i}
$$
and
$$\psi (y)= \sum_{j=1}^{N-1}f(M-1,j)y^{j}= \sum_{j=1}^{N-1}f(1,j)y^{j}.$$
Then,
$$
(4-x-x^{-1}-y-y^{-1})F(x,y) =\\
{\frac {\left (4\,x-4\,x^{M}\right )\left (y-y^{N}\right )}{\left (x-1
\right )\left (y-1\right )}} - (x^{M}+1)\psi (y) - (y^{N}+1)\phi (x). 
$$
So,
$$F(x,y)= $$
$${\frac {xy\left [4\,\left (x-x^{M}\right )\left (y-y^{N}\right )-
\left (x-1\right )\left (y-1\right )\left (x^{M}+1\right )\psi(y)-
\left (x-1\right )\left (y-1\right )\left (y^{N}+1\right )\phi(x)
\right ]}{\left (x-1\right )\left (y-1\right )\left (4xy-x^{2}y-xy^{2}-x
-y\right )}}.$$

Now $\phi (x)$ and $\psi (y)$ are both polynomials, so let
$$\phi (x) = a_{1}x+a_{2}x^{2}+ \dots +a_{M-1}x^{M-1}$$
and
$$\psi (y) = b_{1}y+b_{2}y^{2}+ \dots +b_{N-1}y^{N-1}.$$

We want the numerator of $F(x,y)$ to be $0$ modulo the denominator. To find the
values of $\phi (x)$ and $\psi (y)$ which satisfy this condition, take
the normal form of the numerator of $F(x,y)$ with respect to the ( one
element Gr\" obner basis of the) denominator. Then set each coefficient
of the reduced polynomial to
$0$ and solve the set of equations. From this, $\phi$ and $\psi$,
and hence the generating function, can be obtained.


\section*{Applications to Other Problems}
This procedure has been incorporated in a {\sc maple} program. 
Many other problems involving two-dimensional random walks can be, and
have been, solved by {\sc maple} programs based on the above method. For
example, consider an $[0,M] \times [0,N]$ grid. A particle can move north,
south, east, or west with probabilities $p_{1},p_{2},p_{3}$, and $p_{4}$
respectively. Along the borders there are points where the particle can
escape and points where the particle is captured. We wrote a {\sc maple}
program which automatically computes the preliminary generating function
for the escape probabilities,
then uses Gr\" obner basis functions to find the true generating function. A
program has also been written which deals with more general (not
necessarily convex) shapes. The preliminary manipulations are almost 
completely mechanized. Moreover, this method could also be used for even 
higher dimensions, which would reduce the problem to a smaller system of 
equations ($O(N^{r-1})$ as opposed to $O(N^{r})$ for an $r$-dimensional 
problem).

\section*{Notes}

The programs mentioned above, along with others, may be obtained by
e-mail from {\it craig@euclid.math.temple.edu} upon request.

\begin{thebibliography}{99}

\bibitem{bi : nsee} Birkhoff,G. (1972).The Numerical Solutions of
Elliptic Equations, {\em Regional Conference Series in Applied
Mathematics} {\bf 1}, SIAM: Philadelphia.

\bibitem{bu : tbrpcf} Buchberger, B. (1976).Theoeretical Basis for the
Reduction of Polynomials to Canonical Forms. {\em SIGSAM Bulletin.}
{\bf 40}, 19-29.

\bibitem{dst : ca} Davenport,J.H., Siret, Y., and Tournier, E. (1988).
{\em Computer Algebra: Systems and Algorithms for Algebraic Computation}.
Academic Press: London.

\bibitem{ds : rwaen} Doyle, P.D., Snell, J.L. (1984). {\em Random Walks and
Electrical Networks}. The Mathematical Association of America.

\bibitem{e : ttp} Engle, A. (1993).  The Computer Solves the Three Tower
Problem, {\em American Math. Monthly} {\bf 100}, 62--64.

\bibitem{f : iptia} Feller, W.D. (1950). {\em An Introduction to Probability
Theory and Its Applications}, New York: John Wiley and Sons.

\bibitem{z : mct} Zeilberger, D. (1991). The Method of Creative
Telescoping, {\em J. Symbolic Computation} {\bf 11}, 195--204.

\end{thebibliography}

\end{document}




