Gambler’s Ruin Problem in Several Dimensions

Andrej Kmet and Marko Petkovšek^{∗}

Department of Mathematics and Mechanics Faculty of Mathematics and Physics University of Ljubljana Jadranska 19, SI-1000 Ljubljana, Slovenia

                            Andrej.Kmet@fmf.uni-lj.si
                          Marko.Petkovsek@fmf.uni-lj.si

Abstract

We give an explicit solution of the Gambler’s Ruin Problem when the players use two or more currencies. We also determine the asymptotics of the expected duration of the game when both players have equal amounts of each currency.

1 Introduction

In the one-dimensional Gambler’s Ruin Problem two players start out with i and N − i dollars, respectively. At each step they toss a fair coin to decide who wins a dollar from the opponent. The game is over when one of them goes bankrupt. It is well known that the expected duration of the game is i(N − i) (see [3], or almost any other textbook on probability).

In the two-dimensional variant [6] the players use two different currencies, say dollars and euros. They start out with (i dollars, j euros) and (N − i dollars, M − j euros), respectively. At each step they toss fair coins to decide the currency and the winner. The game is over when one of them runs out of either currency. What is the expected duration of the game? According to [6], no closed-form solution of this problem is known to exist, and probably doesn’t exist.

Denote by game(i, j) the game with the first player’s initial assets equal to (i, j). Assume that 1 ≤ i ≤ N − 1 and 1 ≤ j ≤ M − 1. Then after the first step, game(i, j) turns into one of game(i + 1, j), game(i − 1, j), game(i, j + 1), or game(i, j − 1), each with probability 1/4. It follows that the expected duration a_{i,j} of game(i, j) satisfies the recurrence equation

a_{i,j} = \frac{a_{i+1,j} + a_{i−1,j} + a_{i,j+1} + a_{i,j−1}}{4} + 1 (1 ≤ i ≤ N − 1, 1 ≤ j ≤ M − 1) (1)

∗Supported in part by MZT RS under grant J2-8549.

and the boundary conditions

a_{0,j} = a_{N,j} = a_{i,0} = a_{i,M} = 0 (0 ≤ i ≤ N, 0 ≤ j ≤ M).  (2)

The unknown a_{i,j}, 1 ≤ i ≤ N − 1, 1 ≤ j ≤ M − 1, can be obtained from (1) by straightforward linear algebra. Instead of solving this linear system of (N − 1)(M − 1) equations, Orr and Zeilberger [6] have shown how to obtain the values a_{1,j} = a_{N−1,j} (1 ≤ j ≤ M − 1) and a_{i,1} = a_{i,M−1} (1 ≤ i ≤ N − 1) from a system containing O(N + M) equations only. The remaining a_{i,j} can then be computed directly from the recurrence

a_{i+1,j} = 4 a_{i,j} − (a_{i−1,j} + a_{i,j+1} + a_{i,j−1}) − 4 (1 ≤ i ≤ N − 3, 2 ≤ j ≤ M − 2). (3)

Because of the obvious symmetries a_{i,j} = a_{N−i,j} = a_{i,M−j} = a_{N−i,M−j}, it suffices to compute a quarter of these numbers only. Note however that in floating-point arithmetic, this computation is numerically unstable due to the coefficient 4 in front of a_{i,j} in (3).

In this paper we give an explicit solution of the two-dimensional Gambler’s Ruin Problem (Sec. 2) as well as of its obvious generalization to larger numbers of currencies (Sec. 3). Although our solution is not in closed form because it contains a double sum, it provides a direct way to compute the expected duration of the game without having to solve linear systems or use recursive computations. In Section 4 we express the middle element a_{N/2, N/2} as a single sum when N is a power of two. In Section 5 we determine the asymptotics of the middle element a_{N/2, N/2, ..., N/2} (when N is even) for any number of currencies.

2 Solution of the 2-dimensional problem

For the sake of simplicity we henceforth assume that M = N, but it is straightforward to generalize our results to the case M ≠ N. Let A = [a_{i,j}]_{i,j=1}^{N−1} be the matrix of unknown values a_{i,j}. Writing (1) in the form

(a_{i,j−1} − 2a_{i,j} + a_{i,j+1}) + (a_{i−1,j} − 2a_{i,j} + a_{i+1,j}) = −4 (1 ≤ i, j ≤ N − 1)

and using (2) we see that A satisfies the matrix equation

AD + DA = −4J  (4)

where

D = [ \begin{matrix}−2 & 1 & 0 & 0 & … & 0 \\ 1 & −2 & 1 & 0 & … & 0 \\ 0 & 1 & −2 & 1 & … & 0 \\ ⋮ & ⋮ & ... & ... & ... & ⋮ \\ 0 & 0 & … & 1 & −2 & 1 \\ 0 & 0 & … & 0 & 1 & −2 \end{matrix} ]

is an (N −1)×(N −1) symmetric tridiagonal Toeplitz matrix, and J is the (N −1)×(N −1) matrix of ones. This equation can be solved explicitly in terms of the eigenvalues and eigenvectors of D.

Indeed, let Λ = diag(λ_{k}) be the diagonal matrix of the eigenvalues of D, and P the orthogonal matrix whose columns are the corresponding eigenvectors of D. Then D = P ΛP^{T} and P^{T} P = I. Multiplying (4) with P^{T} on the left and with P on the right gives

BΛ + ΛB = R

where B = P^{T} AP, and

R = −4 P^{T} JP.  (5)

Then

b_{i,j} = \frac{r_{i,j}}{λ_{i} + λ_{j}} (1 ≤ i, j ≤ N − 1),  (6)

so, using (5) and (6), the unknown matrix

A = P BP^{T}  (7)

can be expressed explicitly in terms of Λ and P. These, however, are well known, namely,

λ_{k} = −4 sin^{2} \frac{kπ}{2N} (1 ≤ k ≤ N − 1),  (8)
P_{i,j} = \sqrt{\frac{2}{N}} sin \frac{πij}{N} (1 ≤ i, j ≤ N − 1)  (9)

(cf. [2, Sec. 2.6, item 7^{0}]). Therefore by (5),

r_{i,j} = −4 ∑_{l=1}^{N−1} ∑_{k=1}^{N−1} P_{l,i}P_{k,j} = −\frac{8}{N} ∑_{l=1}^{N−1} sin \frac{πil}{N} ∑_{k=1}^{N−1} sin \frac{πjk}{N} (1 ≤ i, j ≤ N − 1).

This sum is easy to evaluate:

∑_{k=1}^{N−1} sin \frac{πjk}{N} = ℑ ∑_{k=0}^{N−1} e^{\frac{kjπı}{N}} = ℑ \frac{e^{jπı} − 1}{e^{\frac{jπı}{N}} − 1} = \{ \begin{matrix}0, & j even, \\ cot \frac{jπ}{2N} , & j odd, \end{matrix}  (10)

so

r_{i,j} = \{ \begin{matrix}0, & ij even, \\ −\frac{8}{N} cot \frac{iπ}{2N} cot \frac{jπ}{2N} , & ij odd. \end{matrix}

From (6) and (8) it follows that

b_{i,j} = \{ \begin{matrix}0, & ij even, \\ \frac{2 cot \frac{iπ}{2N} cot \frac{jπ}{2N}}{N(sin^{2} \frac{iπ}{2N} +sin^{2} \frac{jπ}{2N})}, & ij odd. \end{matrix}

Finally we obtain from (7) and (9)

a_{i,j} = \frac{4}{N^{2}} ∑_{{k=1 \atop k odd}}^{N−1} sin \frac{jkπ}{N} cot \frac{kπ}{2N} ∑_{{l=1 \atop l odd}}^{N−1} \frac{sin \frac{ilπ}{N} cot \frac{lπ}{2N}}{sin^{2} \frac{kπ}{2N} + sin^{2} \frac{lπ}{2N}} (0 ≤ i, j ≤ N).  (11)

This is the promised explicit expression for the expected duration of game(i, j). It can also be obtained by applying the Discrete Fourier Transform to (1) (see, e.g., [5, Thm. 3a]). Note that it requires only O(N^{2}) time to compute a_{i,j} given i and j.

In the case i = 1, the method of [7] yields the explicit formula

a_{1,j} = \frac{4}{N} ∑_{k=0}^{N−1} sin \frac{(2k + 1)jπ}{N} cot \frac{(2k + 1)π}{2N} \frac{ξ^{N} − ξ}{(ξ^{N} + 1)(ξ − 1)}

where

ξ = (sin \frac{(2k + 1)π}{2N} − \sqrt{1 + sin^{2} \frac{(2k + 1)π}{2N}} )^{2} .

This indicates the possibility that a_{i,j} could be given in the form of a single sum.

3 Generalization to higher dimensions

In d dimensions, the two players use d different currencies. Their initial assets are (i_{1}, i_{2}, …, i_{d}) and (N_{1} − i_{1}, N_{2} − i_{2}, …, N_{d} − i_{d}), respectively. At each step they first decide the currency, then the winner, both uniformly at random. The game is over as soon as one of them runs out of any of the currencies. Let a_{i_{1},i_{2},...,i_{d}} denote the expected duration of the game. Again we assume for simplicity that N_{1} = N_{2} = ⋯ = N_{d} = N.

The solution of the general problem can be expressed explicitly by means of tensor product of matrices. If A is any matrix and B is a matrix of order m × n, their tensor product is defined as

A ⊗ B = [ \begin{matrix}b_{1,1}A & b_{1,2}A & ⋯ & b_{1,n}A \\ b_{2,1}A & b_{2,2}A & ⋯ & b_{2,n}A \\ ⋮ & ⋮ & ... & ⋮ \\ b_{m,1}A & b_{m,2}A & ⋯ & b_{m,n}A \end{matrix} ] .

In particular, if I denotes the identity matrix,

A ⊗ I = [ \begin{matrix}A \\ & A \\ & & . & . & . \\ & & & & & A \end{matrix} ] .

Now consider Sylvester’s equation

AD + CA = R

where C, D, R are known matrices, A is an unknown matrix, A and R are of order p × q, C is of order p × p, and D is of order q × q (cf. [4, Sec. 7.6.3]). This is a system of linear equations

∑_{s=1}^{n} a_{i,s}d_{s,j} + ∑_{t=1}^{n} c_{i,t}a_{t,j} = r_{i,j} (1 ≤ i, j ≤ n)  (12)

for the unknown a_{i,j}. Let a resp. r be the vectors obtained by stacking the columns of A resp. R one above another. If the equations are ordered lexicographically first by j and then by i, we can rewrite (12) as a vector equation

(C ⊗ I_{q} + I_{p} ⊗ D^{T}) a = r  (13)

where I_{n} denotes the identity matrix of order n. Using this on Eqn. (4) which is a special case of Sylvester’s equation with C = D, D = D^{T}, and p = q = N − 1, we have

D_{1}a = −4 r  (14)

where D_{1} = D ⊗ I_{N−1} + I_{N−1} ⊗ D and r is a vector with all components equal to 1. Because D_{1} = (P ⊗ P)(Λ ⊗ I_{N−1} + I_{N−1} ⊗ Λ)(P ⊗ P)^{T} is the spectral decomposition of D_{1}, we obtain a compact form of solution of the two-dimensional problem,

a = −4 (P ⊗ P)(Λ ⊗ I_{N−1} + I_{N−1} ⊗ Λ)^{−1}(P ⊗ P)^{T} r.

In the three-dimensional case, we have the recurrence

(a_{i,j,k−1} − 2a_{i,j,k} + a_{i,j,k+1})
+ (a_{i,j−1,k} − 2a_{i,j,k} + a_{i,j+1,k})  (15)
+ (a_{i−1,j,k} − 2a_{i,j,k} + a_{i+1,j,k}) = −6 (1 ≤ i, j, k ≤ N − 1).

Let A = [a_{.,.,1} | a_{.,.,2} | ⋯ | a_{.,.,N−1}] be the (N − 1)^{2} × (N − 1) matrix whose k-th column a_{.,.,k} is obtained by stacking the columns of the (N − 1) × (N − 1) matrix [a_{i,j,k}]_{i,j=1}^{N−1} one above another. Then (15) is equivalent to Sylvester’s equation

AD + D_{1}A = −6J.

Let a be the vector obtained by stacking the columns of A one above another. Using (13) we have

D_{2} a = −6 r,

where D_{2} = D_{1} ⊗ I_{N−1} + I_{(N−1)^{2}} ⊗ D, and r is a vector of ones.

In the d-dimensional case we recursively define three sequences of matrices by setting D_{0} = D, Λ_{0} = Λ, P_{0} = P, and, for k = 1, 2, …, d − 1,

D_{k} = D_{k−1} ⊗ I_{N−1} + I_{(N−1)^{k}} ⊗ D,
Λ_{k} = Λ_{k−1} ⊗ I_{N−1} + I_{(N−1)^{k}} ⊗ Λ,
P_{k} = P_{k−1} ⊗ P.

Then for all k we have P_{k}^{T} P_{k} = I_{(N−1)^{k+1}}, D_{k} = P_{k}Λ_{k}P_{k}^{T}, and

D_{d−1} a = −2d r

where a is the vector with components a_{i_{1},i_{2},...,i_{d}}, ordered lexicographically first by i_{d}, then by i_{d−1}, …, last by i_{1}, and r is a vector of ones. Thus the solution in the d-dimensional case is given by

a = −2d P_{d−1}Λ_{d−1}^{−1}P_{d−1}^{T} r.  (16)

4 A faster way to solution for N = 2^{p}

Let d = 2 again, and N = 2^{p} for some integer p > 0. In this case the system (14) can be solved efficiently using the block cyclic reduction technique or Buneman’s algorithm (see [1], [8], or [4, Sec. 4.5.4]). Recall that D_{1} = D⊗I_{N−1}+I_{N−1}⊗D, so (14) has the (N −1)×(N −1) block-tridiagonal form

[ \begin{matrix}G & I & 0 & 0 & … & 0 \\ I & G & I & 0 & … & 0 \\ 0 & I & G & I & … & 0 \\ ⋮ & ⋮ & ... & ... & ... & ⋮ \\ 0 & 0 & … & I & G & I \\ 0 & 0 & … & 0 & I & G \end{matrix} ] [ a_{1} ; a_{2} ; a_{3} ; ⋮ ; a_{N−2} ; a_{N−1} ] = [ b ; b ; b ; ⋮ ; b ; b ]

where G = D − 2I is of size (N − 1) × (N − 1), a_{j} is the j-th column of A, and b is a vector with all components equal to −4.

Multiplying every other row by −G and adding the two neighboring rows to it we eliminate the odd-numbered blocks and obtain the (N/2 − 1) × (N/2 − 1) block-tridiagonal system

[ \begin{matrix}G_{2} & I & 0 & 0 & … & 0 \\ I & G_{2} & I & 0 & … & 0 \\ 0 & I & G_{2} & I & … & 0 \\ ⋮ & ⋮ & ... & ... & ... & ⋮ \\ 0 & 0 & … & I & G_{2} & I \\ 0 & 0 & … & 0 & I & G_{2} \end{matrix} ] [ a_{2} ; a_{4} ; a_{6} ; ⋮ ; a_{N−4} ; a_{N−2} ] = [ b_{2} ; b_{2} ; b_{2} ; ⋮ ; b_{2} ; b_{2} ]

where G_{2} = 2I − G^{2} and b_{2} = (2I − G)b. Repeating this procedure p − 1 times we come to a single-block system

G_{p}a_{N/2} = b_{p}  (17)

where we recursively defined G_{1} = G, b_{1} = b, G_{k+1} = 2I − G_{k}^{2}, and b_{k+1} = (2I − G_{k})b_{k}. Let

Λ_{1} = Λ − 2I,
Λ_{k+1} = 2I − Λ_{k}^{2} for k ≥ 0  (18)

where D = P ΛP^{T}, P^{T} P = I, and the elements of Λ and P are given in (8) resp. (9). By induction on k it follows that G_{k} = P Λ_{k}P^{T} and b_{k} = P (2I + Λ_{k})Λ^{−1}P^{T} b. Hence from (17)

a_{N/2} = P (I + 2Λ_{p}^{−1})Λ^{−1}P^{T} b.  (19)

The remaining a_{j} can now be found by back-substitution from the intermediate systems. For the middle element a_{N/2,N/2} = a_{2^{p−1}, 2^{p−1}} we obtain from (19), (8), (9), and (10)

a_{N/2,N/2} = \frac{2}{N} ∑_{k=0}^{N/2−1}(−1)^{k} (1 + \frac{2}{λ_{2k+1}^{(p)}} ) \begin{matrix}cos & \frac{(2k+1)π}{2N} \\ sin & 3 & \frac{(2k+1)π}{2N} \end{matrix}  (20)

where λ_{j}^{(k)}, the j-th element of Λ_{k}, satisfies

λ_{j}^{(1)} = λ_{j} − 2,
λ_{j}^{(k+1)} = 2 − (λ_{j}^{(k)})^{2} for k ≥ 0.

Using the substitution λ_{j}^{(k)} = −2 cosh a_{k} we find a_{k+1} = 2a_{k} and a_{1} = cosh^{−1}(1 − λ_{j}/2), so a_{k} = 2^{k−1}a_{1} and λ_{j}^{(k)} = −2 cosh (2^{k−1} cosh^{−1}(1 − λ_{j}/2)). With (8) and (20) we finally obtain

a_{N/2,N/2} = \frac{2}{N} ∑_{k=0}^{N/2−1}(−1)^{k} (1 − \frac{1}{cosh (\frac{N}{2} cosh^{−1} (1 + 2 sin^{2} \frac{(2k+1)π}{2N} ))}) \begin{matrix}cos & \frac{(2k+1)π}{2N} \\ sin & 3 & \frac{(2k+1)π}{2N} \end{matrix} .

Note that this formula requires only O(N) operations to compute a_{N/2,N/2}.

5 Some asymptotics

Clearly the expected duration of the game is largest when the initial assets of the two players in each currency are as close to each other as possible. Here we assume that N_{1} = N_{2} = ⋯ = 2n, and find an asymptotic formula for the middle element a_{n,n,...,n}.

Theorem 1 Let d = 2. If N_{1} = N_{2} = 2n then

a_{n,n} = c_{2} n^{2} + O(log^{2} n)

where

c_{2} = \frac{256}{π^{4}} ∑_{k=0}^{∞} ∑_{l=0}^{∞} \frac{(−1)^{k+l}}{(2k + 1)(2l + 1) ((2k + 1)^{2} + (2l + 1)^{2})} ≈ 1.17874.  (21)

Proof: From (11) we have

a_{n,n} = \frac{1}{n^{2}} ∑_{{k=1 \atop k odd}}^{2n−1} ∑_{{l=1 \atop l odd}}^{2n−1} \frac{sin \frac{kπ}{2} sin \frac{lπ}{2} cot \frac{kπ}{4n} cot \frac{lπ}{4n}}{sin^{2} \frac{kπ}{4n} + sin^{2} \frac{lπ}{4n}}

= \frac{1}{n^{2}} ∑_{k=0}^{n−1} ∑_{l=0}^{n−1} \frac{(−1)^{k+l} cot \frac{(2k+1)π}{4n} cot \frac{(2l+1)π}{4n}}{sin^{2} \frac{(2k+1)π}{4n} + sin^{2} \frac{(2l+1)π}{4n}} .  (22)

Denote by s_{n} the expression obtained by replacing sin x with x and cos x with 1 everywhere in (22):

s_{n} = \frac{256}{π^{4}}n^{2} ∑_{k=0}^{n−1} ∑_{l=0}^{n−1} \frac{(−1)^{k+l}}{(2k + 1)(2l + 1) ((2k + 1)^{2} + (2l + 1)^{2})}.  (23)

For 0 < x, y < π/2 we have

\frac{(1 − \frac{x^{2}}{2} ) (1 − \frac{y^{2}}{2} )}{xy(x^{2} + y^{2})} < \frac{cot x cot y}{sin^{2} x + sin^{2} y} < \frac{1}{xy(x^{2} + y^{2})}.

The first inequality follows from cos x > 1−x^{2}/2 and sin x < x, the second from the fact that the Taylor coefficients of f(x, y) = (sin^{2} x + sin^{2} y) tan x tan y − xy(x^{2} + y^{2}) are nonnegative at x = y = 0. Thus

0 < \frac{1}{xy(x^{2} + y^{2})} − \frac{cot x cot y}{sin^{2} x + sin^{2} y} < \frac{1 − (1 − \frac{x^{2}}{2} ) (1 − \frac{y^{2}}{2} )}{xy(x^{2} + y^{2})} = \frac{x^{2} + y^{2} − \frac{x^{2}y^{2}}{2}}{2xy(x^{2} + y^{2})} < \frac{1}{2xy}

which shows that the error committed in replacing \frac{cot x cot y}{sin^{2} x+sin^{2} y} with \frac{1}{xy(x^{2}+y^{2})} is bounded by \frac{1}{2xy}. It follows that

|a_{n,n} − s_{n}| < \frac{1}{n^{2}} ∑_{k=0}^{n−1} ∑_{l=0}^{n−1} \frac{1}{2\frac{(2k+1)π}{4n} \frac{(2l+1)π}{4n}} = \frac{8}{π^{2}} ∑_{k=0}^{n−1} \frac{1}{2k + 1} ∑_{l=0}^{n−1} \frac{1}{2l + 1} = \frac{8}{π^{2}} (H_{2n} − \frac{1}{2}H_{n})^{2}

where H_{n} = ∑_{k=1}^{n} 1/k is the n-th harmonic number. As ln n + γ < H_{n} ≤ ln n + 1, it follows that a_{n,n} − s_{n} = O(log^{2} n).

Next replace the sum in (23) by the infinite sum

t_{n} = \frac{256}{π^{4}}n^{2} ∑_{k=0}^{∞} ∑_{l=0}^{∞}(−1)^{k+l}r_{k,l} = c_{2} n^{2}  (24)

where

r_{k,l} = \frac{1}{(2k + 1)(2l + 1) ((2k + 1)^{2} + (2l + 1)^{2})}.

The sum in (24) is absolutely convergent, therefore it can be rearranged into

t_{n} = \frac{256}{π^{4}}n^{2} ∑_{j=0}^{∞}(−1)^{j}d_{j}  (25)

where

d_{j} = ∑_{k=0}^{j} r_{k,j−k}

is the absolute value of the sum of the j-th diagonal of the original array. We transform s_{n} into t_{n} in two steps: first we approximate s_{n} by

u_{n} = \frac{256}{π^{4}}n^{2} ∑_{j=0}^{n}(−1)^{j}d_{j},

then we replace u_{n} by t_{n}.

It is possible to show that d_{j} > d_{j+1} for all j ≥ 0. Therefore the error in truncating the alternating sum (24) after the n-th diagonal does not exceed

\frac{256}{π^{4}}n^{2}d_{n+1} = \frac{256}{π^{4}}n^{2} ∑_{k=0}^{n+1} r_{k,n+1−k} ≤ \frac{256}{π^{4}}n^{2}(n + 2)r_{0,n+1} = \frac{256 n^{2}(n + 2)}{π^{4}(2n + 3)(4n^{2} + 12n + 10)} < \frac{32}{π^{4}}.

Thus both differences s_{n} − u_{n} and u_{n} − t_{n} are O(1), and the theorem is proved. □

Theorem 2 If N_{1} = N_{2} = ⋯ = N_{d} = N = 2n then

a_{n,n,...,n} ∼ c_{d} n^{2}

where

c_{d} = d (1 − \frac{2^{2d+1}}{π^{d+1}} ∑_{k_{1},k_{2},...,k_{d−1}≥0} \frac{(−1)^{∑_{j=1}^{d−1} k_{j}} ∏_{j=1}^{d−1} \frac{1}{2k_{j}+1} ∑_{j=1}^{d−1} \frac{1}{(2k_{j}+1)^{2}}}{cosh \frac{π}{2} \sqrt{∑_{j=1}^{d−1}(2k_{j} + 1)^{2}}} ) .  (26)

Proof: Consider the following Poisson partial differential equation in d dimensions

Δu = −2d  (27)

on the d-dimensional hypercube Q_{d} = [0, 1]^{d}, with zero boundary conditions: u|_{∂Q_{d}} = 0. Following the finite-difference method to obtain an approximate solution, we discretize the problem by introducing the d-dimensional sequence of values on the square grid

a_{i_{1},i_{2},...,i_{d}} = \frac{1}{h^{2}} u (\frac{i_{1}}{N} , \frac{i_{2}}{N} , …, \frac{i_{d}}{N} )

where h = 1/N and 0 ≤ i_{1}, i_{2}, …, i_{d} ≤ N. Replacing second derivatives by second differences divided by h^{2}, Eqn. (27) turns into the d-dimensional Gambler’s Ruin recurrence for a_{i_{1},i_{2},...,i_{d}}. Thus, for large n,

a_{n,n,...,n} ∼ 4n^{2} u (\frac{1}{2}, \frac{1}{2}, …, \frac{1}{2}) .  (28)

It remains to solve Eqn. (27) exactly. First we homogenize the equation by introducing

v(x_{1}, x_{2}, …, x_{d}) = u(x_{1}, x_{2}, …, x_{d}) + ∑_{i=1}^{d} x_{i}(x_{i} − 1)  (29)

which satisfies

Δv = 0, v|_{x_{i}=0} = v|_{x_{i}=1} = ∑_{{j=1 \atop j≠i}}^{d} x_{j}(x_{j} − 1).

Then we solve this Dirichlet problem exactly by the standard Fourier-series method and obtain

v(x_{1}, x_{2}, …, x_{d}) = −\frac{2^{2d−1}}{π^{d+1}} ∑_{{k_{1},...,k_{d−1}≥0 \atop k_{1},...,k_{d−1} odd}} \begin{matrix}∏ & {d−1 \atop j=1} & \frac{1}{k_{j}} & ∑ & {d−1 \atop j=1} & \frac{1}{k_{j}^{2}} \\ cosh & \frac{π}{2} & \sqrt{∑_{j=1}^{d−1} k_{j}^{2}} \end{matrix} ∑_{i=1}^{d−1} f_{i}(x_{1}, …, x_{d}; k_{1}, …, k_{d−1}),

(30)

where

f_{i}(x_{1}, …, x_{d}; k_{1}, …, k_{d−1}) = ∏_{j=1}^{i−1} sin πk_{j}x_{j} ∏_{j=i}^{d−1} sin πk_{j}x_{j+1} cosh (π (x_{i} − \frac{1}{2}) \sqrt{∑_{j=1}^{d−1} k_{j}^{2}}) .

Combining (28) with (29) and (30) yields (26). □

In particular, for d = 2 we obtain from (26) a single-sum expression for c_{2},

c_{2} = 2 (1 − \frac{32}{π^{3}} ∑_{k≥0} \frac{(−1)^{k}}{(2k + 1)^{3} cosh \frac{π}{2}(2k + 1)})

(cf. Eqn. (21)). Some values of c_{d} are shown in Table 1.

\begin{matrix}d & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ c_{d} & 1 & 1.17874 & 1.34911 & 1.51271 & 1.67071 & 1.82399 & 1.97321 & 2.11891 & 2.26149 \end{matrix}

Table 1: Approximate values of c_{d}

References

[1] B. L. Buzbee, G. H. Golub, C. W. Nielson, On direct methods for solving Poisson’s equation, SIAM J. Numer. Anal. 7 (1970) 627–656.

[2] D. M. Cvetković, M. Doob, H. Sachs, Spectra of Graphs, 3rd ed., Johann Ambrosius Barth Verlag, Heidelberg – Leipzig 1995.

[3] W. D. Feller, An Introduction to Probability Theory and its Applications, John Wiley & Sons, New York 1950.

[4] G. H. Golub and C. F. van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University Press, Baltimore 1996.

[5] P. Henrici, Poisson’s equation in a hypercube: Discrete Fourier methods, eigenfunction expansions, Padé approximation to eigenvalues, in: G. H. Golub, Ed., Studies in Numerical Analysis, vol. 24 of MAA Studies in Mathematics, The Mathematical Association of America, 1984, 371–411.

[6] C. R. Orr and D. Zeilberger, A computer algebra approach to the discrete Dirichlet problem, J. Symb. Comp. 18 (1994) 87–90.

[7] M. Petkovšek, Discrete boundary-value problems, in preparation.

[8] P. N. Swarztrauber, Fast Poisson solvers, in: G. H. Golub, Ed., Studies in Numerical Analysis, vol. 24 of MAA Studies in Mathematics, The Mathematical Association of America, 1984, 319–370.