Electronic Journal of Qualitative Theory of Differential Equations Spec. Ed. I, 2009 No. 4, 1–12; http://www.math.u-szeged.hu/ejqtde/

GREEN’S FUNCTION OF A CENTERED PARTIAL DIFFERENCE EQUATION

Richard I. Avery ^{1} and Douglas R. Anderson ^{2}

1 College of Arts and Sciences, Dakota State University Madison, SD 57042 USA e-mail: rich.avery@dsu.edu 2 Department of Mathematics and Computer Science, Concordia College Moorhead, MN 56562 USA e-mail: andersod@cord.edu

Honoring the Career of John Graef on the Occasion of His Sixty-Seventh Birthday

Abstract

Applying a variation of Jacobi iteration we obtain the Green’s function for the centered partial difference equation

Δ_{ww}u(x_{w−1}, y_{z}) + Δ_{zz}u(x_{w}, y_{z−1}) + f(u(x_{w}, y_{z})) = 0,

which is the result of applying the finite difference method to an associated nonlinear partial differential equation of the form

u_{xx} + u_{yy} + h(u) = 0.

We show that approximations of the partial differential equation can be found by applying fixed point theory instead of the standard techniques associated with solving a system of nonlinear equations.

Key words and phrases: Finite difference method, partial difference equations, Green’s function. AMS (MOS) Subject Classifications: 39A12, 37H10.

1 Introduction

The finite difference method is a classical technique used to approximate the solution of a partial differential equation. An application of the technique yields a nonlinear system of equations which can be considered as a nonlinear partial difference equation. Instead of attempting to solve or approximate the solution of the nonlinear system, we will convert the partial difference equation to a fixed point problem whose inversion involves the Green’s function of the partial difference equation. The Green’s function can be interpreted combinatorially through an adjacency matrix. Finding the Green’s function for a partial difference equation has been done; for example, see [1, 2, 3, 4, 5, 6, 7, 8, 12, 13], however none of the techniques are simple nor are the formulas easy to work with when finding solutions which correspond to approximations to partial differential equations. Also, it is difficult to find bounds on the Green’s function to establish the existence of solutions for the classical partial difference equation when these techniques are used. In this paper we apply an elementary iterative technique based on Jacobi iteration to obtain the Green’s functions. A technique that is easy to computationally apply to approximate the solution of a partial differential equation using partial difference equations techniques. Since the Green’s function involves the sum of nonnegative terms (with a combinatorial interpretation) it is easy to establish elementary bounds which can be used in existence of solutions arguments. A basic understanding of difference equations, Green’s functions, and discretization of partial differential equations is assumed; see Kelley and Peterson [10] for background results.

2 Preliminaries

For designation purposes, we will employ interval notation to denote sets of integers, such as [1, N] = {1, 2, …, N}, etc.

The partial difference equation

Δ_{ww}u(x_{w−1}, y_{z}) + Δ_{zz}u(x_{w}, y_{z−1}) + f(u(x_{w}, y_{z})) = 0,

when w ∈ [1, N − 1] and z ∈ [1, N − 1], with boundary conditions

u(x, y) = 0 for all x = x_{0}, x = x_{N}, y = y_{0}, or y = y_{N},

is the result of applying the finite element method to the nonlinear partial differential equation

u_{xx} + u_{yy} + h(u) = 0,

with boundary conditions u(0, y) = 0, 0 ≤ y ≤ 1, and u(x, 0) = 0, 0 ≤ x ≤ 1, where N^{2}f(u) = h(u), x_{w} = \frac{w}{N}, y_{z} = \frac{z}{N}. Thus, this is a problem that is of interest to a large audience extending beyond mathematicians. Note, the arguments throughout can be modified to suit any bounded domain. The Green’s function is the function whose domain is [0, N]^{2} × [1, N − 1]^{2} such that for all (r, s) ∈ [1, N − 1]^{2},

g(w, z, r, s) = 0, if (w, z) ∈ ∂[0, N]^{2},

and

Δ_{ww}g(w − 1, z, r, s) + Δ_{zz}g(w, z − 1, r, s) = −δ_{w,r}δ_{z,s}

for all w ∈ [1, N − 1] and z ∈ [1, n − 1]. Applying classical techniques it can be shown that any solution of the partial difference equation is a fixed point of the operator T defined by

T u(x_{w}, y_{z}) = ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)f(u(x_{r}, y_{s})).

3 Finding the Green’s Function

For each fixed (r, s) ∈ [1, N − 1]^{2}, the Green’s function is a solution of

Δ_{ww}g(w − 1, z, r, s) + Δ_{zz}g(w, z − 1, r, s) = −δ_{w,r}δ_{z,s},

when w ∈ [1, N − 1] and z ∈ [1, N − 1], which we can write as

g(w, z, r, s) = (\frac{1}{4}) [g(w − 1, z, r, s) + g(w + 1, z, r, s)
+ g(w, z − 1, r, s) + g(w, z + 1, r, s) + δ_{w,r}δ_{z,s}].

Thus, for each fixed (r, s) ∈ [1, N − 1]^{2}, the Green’s function is a solution of a linear system of equations which we can write as

G_{r,s} = (\frac{A}{4}) G_{r,s} + \frac{D_{r,s}}{4} .

The column matrix G_{r,s} consists of the Green’s function terms, that is, for each (r, s) ∈ [1, N −1]^{2}, the element in the ith row of G_{r,s} is g(w, z, r, s), where i = (w−1)(N −1)+z. Similarly, for each (r, s) ∈ [1, N − 1]^{2}, the element in the ith row of D_{r,s} is δ_{w,r}δ_{z,s}, where i = (w − 1)(N − 1) + z. The matrix A is the (N − 1)^{2} × (N − 1)^{2} adjacency matrix for the Laplacian molecule. That is, if

i = (w_{i} − 1)(N − 1) + z_{i} and j = (w_{j} − 1)(N − 1) + z_{j},

then the (i, j)th entry of A is 1, if (w_{i}, z_{i}) is adjacent to (w_{j}, z_{j}), and 0 otherwise. The vertices (w_{i}, z_{i}) and (w_{j}, z_{j}) are adjacent if

w_{i} = w_{j} and |z_{i} − z_{j}| = 1, or z_{i} = z_{j} and |w_{i} − w_{j}| = 1.

Since A is an adjacency matrix, for all natural numbers k, the (i, j)th entry of A^{k} is the number of distinct paths of length k from node i to node j that never leave the bounded lattice [1, N − 1]^{2} (see [11] for details), where node i corresponds to the point (w_{i}, z_{i}) in the lattice, and the node j corresponds to the point (w_{j}, z_{j}) in the lattice, again where

i = (w_{i} − 1)(N − 1) + z_{i} and j = (w_{j} − 1)(N − 1) + z_{j}.

For each natural number t and (i, j) ∈ [1, (N − 1)^{2}]^{2}, let a_{i,j}^{(t)} represent the i, jth entry of A^{t}. Then the maximum row sum matrix norm is defined by

∥A^{t}∥_{∞} = max_{i∈[1,(N−1)^{2}]} ∑_{j=1}^{(N−1)^{2}} |a_{i,j}^{(t)}|.

In the theorem below we verify that I − \frac{A}{4} is invertible.

Theorem 3.1 If A is the adjacency matrix corresponding to the centered Laplacian with zero boundary conditions, then I − \frac{A}{4} is invertible. Moreover,

(I − \frac{A}{4})^{−1} = ∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} .

Proof. Let i ∈ [1, (N − 1)^{2}] and (w_{i}, z_{i}) be the lattice point corresponding to the i^{th} node. For any natural number t, the number of walks in an infinite lattice of length t from node i is

4^{t} = (1 + 1 + 1 + 1)^{t} = ∑_{r=0}^{t} ∑_{l=0}^{t−r} ∑_{u=0}^{t−r−l} ( \begin{matrix}t \\ r \end{matrix} ) ( \begin{matrix}t − r \\ l \end{matrix} ) ( \begin{matrix}t − r − l \\ u \end{matrix} ) ,

and if t ≥ ⌈\frac{N}{2} ⌉, then at least one of these walks has left the bounded lattice [1, N − 1]^{2}, since a walk in the direction (all the steps either left, right, up or down) of the closest boundary will leave the bounded lattice [1, N − 1]^{2}. Hence, if t ≥ ⌈\frac{N}{2} ⌉, then for all i ∈ [1, N − 1]^{2},

∑_{j=1}^{(N−1)^{2}} |a_{i,j}^{(t)}| ≤ 4^{t} − 1;

that is, the number of walks from node i to any other node in the bounded lattice in t steps that does not leave the bounded lattice is at most 4^{t} − 1. Therefore,

∥\frac{A^{⌈ \frac{N}{2} ⌉}}{4^{⌈ \frac{N}{2} ⌉}} ∥_{∞} ≤ \frac{4^{⌈ \frac{N}{2} ⌉} − 1}{4^{⌈ \frac{N}{2} ⌉}}

which implies that

I − \frac{A^{⌈ \frac{N}{2} ⌉}}{4^{⌈ \frac{N}{2} ⌉}}

is invertible. Since

I − \frac{A^{⌈ \frac{N}{2} ⌉}}{4^{⌈ \frac{N}{2} ⌉}} = (I − \frac{A}{4}) (I + \frac{A}{4} + \frac{A^{2}}{4^{2}} + ⋯ + \frac{A^{⌈ \frac{N}{2} ⌉−1}}{4^{⌈ \frac{N}{2} ⌉−1}} )

the invertibility of I − \frac{A^{⌈ \frac{N}{2} ⌉}}{4^{⌈ \frac{N}{2} ⌉}} guarantees the invertibility of I − \frac{A}{4} as well as the invert-

ibility of I + \frac{A}{4} + \frac{A^{2}}{4^{2}} + ⋯ + \frac{A^{⌈ \frac{N}{2} ⌉−1}}{4^{⌈ \frac{N}{2} ⌉−1}} . Moreover, it is well known that if (I − \frac{A}{4})^{−1} exists, then

(I − \frac{A}{4})^{−1} = ∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} .

Now that we have that I − \frac{A}{4} is invertible and that

(I − \frac{A}{4})^{−1} = ∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} ,

then, for each fixed (r, s) ∈ [1, N − 1]^{2}, the Green’s function is given by

G_{r,s} = (∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} ) \frac{D_{r,s}}{4} .

Note, that given that I − \frac{A}{4} is invertible, then the Green’s function could be found by Jacobi Iteration. That is, if we let X_{0} = \frac{D_{r,s}}{4} , and for n ≥ 1, let

X_{n} = (\frac{A}{4}) X_{n−1} + \frac{D_{r,s}}{4}

then by iteration, we have

X_{n} = (∑_{k=0}^{n} \frac{A^{k}}{4^{k}} ) (\frac{D_{r,s}}{4} )

and hence

G_{r,s} = lim_{n→∞} X_{n}.

Example 3.1 Suppose N = 4. Then, the adjacency matrix A corresponding to the centered Laplacian with zero boundary conditions (4 × 4 lattice with the boundary removed which results in a 3 × 3 lattice) is

A = [ \begin{matrix}0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \end{matrix} ]

and since

G_{r,s} = (∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} ) (\frac{D_{r,s}}{4} ) ,

lower and upper bounds on ∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} will be crucial in existence of solutions arguments as well as iterative arguments to find an approximation of the original partial differential equation problem.

4 Alternative Inversion Technique

In the literature concerning the existence of positive solutions of various boundary value problems, the most common procedure is to apply fixed point theorems to operators analogous to the operator T above. Rather than taking that approach in this paper, we will establish the existence of a fixed point for a related operator S.

Let P be a space of functions defined on [0, N]^{2}, and define

S : P → P

by

Sv(w, z) = f (∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s)) ,

for (w, z) ∈ [0, N]^{2}. We now show that if v ∈ P is a fixed point of S, then

u(w, z) := ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s)

is a fixed point of T (and hence is a solution of the partial difference equation). To see this, assume v ∈ P is a fixed point of S. Then,

T u(w, z) = ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)f(u(r, s))

= ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)f (∑_{τ=1}^{N−1} ∑_{ψ=1}^{N−1} g(r, s, τ, ψ)v(τ, ψ))

= ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)S(v(r, s))

= ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s)
= u(w, z).

Conversely, assume u ∈ P is a fixed point of T. Then,

v(w, z) := f(u(w, z))

is a fixed point of S. To see this, assume u ∈ P is a fixed point of T. Then,

Sv(w, z) = f (∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s))

= f (∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)f(u(r, s)))
= f (T u(w, z))
= f(u(w, z)) = v(w, z).

It follows easily that there is a one-to-one correspondence between fixed points of T and the fixed points of S. By applying the definition, in an appropriate function space context, it can be shown that, if f is continuous, then both S and T are completely continuous.

5 Application

In this section we will state the definitions that are used in the remainder of the paper.

Definition 5.1 Let E be a real Banach space. A nonempty closed convex set P ⊂ E is called a cone if it satisfies the following two conditions:

(i) x ∈ P, λ ≥ 0 implies λx ∈ P;

(ii) x ∈ P, −x ∈ P implies x = 0.

Every cone P ⊂ E induces an ordering, ≤, in E given by

x ≤ y if and only if y − x ∈ P.

Definition 5.2 An operator is called completely continuous if it is continuous and maps bounded sets into pre-compact sets.

Definition 5.3 Let P be a cone in a real Banach space E and D ⊆ E. Then the operator A : D → E is said to be increasing on D provided x_{1}, x_{2} ∈ D with x_{1} ≤ x_{2} implies Ax_{1} ≤ Ax_{2}.

Definition 5.4 A cone P of a real Banach space E is said to be normal if there exists a positive constant δ such that ∥x + y∥ ≥ δ for all x, y ∈ P with ∥x∥ = ∥y∥ = 1.

The following theorem is an elementary fact about normal cones. A proof can be found in [9].

Theorem 5.1 Let P be a cone in a real Banach space E. The cone P is normal if and only if the norm of the Banach space E is semi-monotone; that is, there exists a constant N > 0 such that 0 ≤ x ≤ y implies that ∥x∥ ≤ N∥y∥.

The next two theorems concerning the convergence of Picard iterates are restatements of theorems that can be found in [9].

Theorem 5.2 Let P be a normal cone in a real Banach space E and A : P → E be a completely continuous operator. If u ∈ P with Au ≤ u and there exists a v ∈ P such that v ≤ A^{n}u for all n ∈ N and A is increasing on [v, u], then

{A^{n}u}_{n=1}^{∞}

is a decreasing sequence bounded below by v ∈ P, and there exists a fixed point u^{∗} ∈ P of A such that

u^{∗} = lim_{n→∞} A^{n}u

with

v ≤ u^{∗} ≤ A^{n}u ≤ A^{n−1}u ≤ ⋯ ≤ Au.

Theorem 5.3 Let P be a normal cone in a real Banach space E and A : P → E be a completely continuous operator. If u ∈ P with u ≤ Au and there exists a v ∈ P such that v ≥ A^{n}u for all n ∈ N and A is increasing on [u, v], then

{A^{n}u}_{n=1}^{∞}

is an increasing sequence bounded above by v ∈ P, and there exists a fixed point u^{∗} ∈ P of A such that

u^{∗} = lim_{n→∞} A^{n}u

with

v ≥ u^{∗} ≥ A^{n}u ≥ A^{n−1}u ≥ ⋯ ≥ Au.

We now present our solutions result as a fixed point application.

Theorem 5.4 Suppose there exist positive real numbers m and M, with 0 < m < M, and an increasing continuous function f : [m, ∞) → [0, ∞) such that

f(x) ≤ M for x ≤ max_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s) M

and

m ≤ f(x) for min_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)m ≤ x.

Then,

v_{l} = lim_{n→∞} v_{l,n} and v_{u} = lim_{n→∞} v_{u,n}

are fixed points of S, where for all natural numbers n, v_{l,n+1} = S(v_{l,n}) and v_{u,n+1} = S(v_{u,n}), with

v_{l,0}(w, z) = \{ \begin{matrix}m, & (w, z) ∈ [1, N −1]^{2}, \\ f(0), & otherwise, \end{matrix}

and

v_{u,0}(w, z) = \{ \begin{matrix}M, & (w, z) ∈ [1, N −1]^{2}, \\ f(0), & otherwise. \end{matrix}

Moreover, this implies that

u^{∗}(w, z) = ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{u}(r, s)

and

u_{∗}(w, z) = ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{l}(r, s)

are fixed points of the operator T and are therefore solutions (maybe equal) of the original partial difference equation.

Proof. Let

E := {u : u : [0, N]^{2} → ℝ}

and

P := {u : u ∈ E and u ≥ 0}.

Then P is a normal cone in the Banach space E with the maximum norm

∥u∥ = max_{(w,z)∈[0,N]^{2}} u(w, z).

For any v ∈ P, with min_{(w,z)∈[1,N−1]^{2}} v(w, z) ≥ m,

min_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)m ≤ min_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s),

and for any v ∈ P, with max_{(w,z)∈[1,N−1]^{2}} v(w, z) ≤ M,

max_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v(r, s) ≤ max_{(w,z)∈[1,N−1]^{2}} ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)M.

Therefore, for all (w, z) ∈ [1, N − 1]^{2}, we have

v_{l,1}(w, z) = S(v_{l,0}(w, z)) ≥ m = v_{l,0}(w, z)

and

v_{u,1}(w, z) = S(v_{u,0}(w, z)) ≤ M = v_{u,0}(w, z).

Note that, if (w, z) ∉ [1, N − 1]^{2} (that is, (w, z) is on the boundary), then

v_{l,1}(w, z) = f(0) = S(v_{l,0}(w, z)) and v_{u,1}(w, z) = f(0) = S(v_{u,0}(w, z)).

Also, for any v_{1}, v_{2} ∈ P with v_{l,0} ≤ v_{1} ≤ v_{2} ≤ v_{u,0}, and for all (w, z) ∈ [0, N]^{2}, we have

∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{1}(r, s) ≤ ∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{2}(r, s).

Thus, since f is increasing,

f (∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{1}(r, s)) ≤ f (∑_{r=1}^{N−1} ∑_{s=1}^{N−1} g(w, z, r, s)v_{2}(r, s)) ,

and hence S(v_{1}) ≤ S(v_{2}). Therefore, by Theorems 5.2 and 5.3, we have that

v_{l} = lim_{n→∞} v_{l,n} and v_{u} = lim_{n→∞} v_{u,n}

are fixed points of S. Moreover, we have that

v_{l,0} ≤ v_{l,1} ≤ ⋯ ≤ v_{l,n} ≤ ⋯ ≤ v_{l} ≤ v_{u} ≤ ⋯ ≤ v_{u,n} ≤ ⋯ ≤ v_{u,1} ≤ v_{u,0}.

Example 5.1 Consider the partial difference equation

Δ_{ww}u (\frac{w − 1}{4}, \frac{z}{4}) + Δ_{zz}u (\frac{w}{4}, \frac{z − 1}{4}) + 10^{2} arctan (u (\frac{w}{4}, \frac{z}{4})) = 0,

when w, z ∈ {1, 2, 3}, with boundary conditions

u(0, t) = u(1, t) = u(t, 0) = u(t, 1) = 0,

for all t = \frac{k}{4}, where k ∈ \{0, 1, 2, 3, 4\}, which is the result of applying the finite difference method to the boundary value problem

u_{xx} + u_{yy} + 40^{2} arctan(u) = 0

on the unit square with zero boundary conditions taking horizontal and vertical steps of size one fourth. For all (w, z) ∈ [1, 3]^{2}, there are at least two adjacent vertices; hence,

∑_{r=1}^{3} ∑_{s=1}^{3} g(w, z, r, s)(40π) ≥ (\frac{1}{4}) (1 + \frac{2}{4}) (40π) = 15π,

which can also be obtained by summing the entries in the ith row of

(∑_{k=0}^{1} \frac{A^{k}}{4^{k}} ) (\frac{1}{4}) ,

where i = w(4 − 1) + z, since as noted in Example 3.1,

G_{r,s} = (∑_{k=0}^{∞} \frac{A^{k}}{4^{k}} ) (\frac{D_{r,s}}{4} ) .

One can find a solution by iteration to the partial difference equation that approximates the solution of the partial differential equation applying Theorem 5.4, with M = 50π and m = 40π, since for all (w, z) ∈ [1, 3]^{2},

40π < 100 arctan(15π) < 100 arctan (∑_{r=1}^{3} ∑_{s=1}^{3} g(w, z, r, s)(40π))

and

100 arctan (∑_{r=1}^{3} ∑_{s=1}^{3} g(w, z, r, s)(50π)) < 100 (\frac{π}{2} ) = 50π.

Using Maple 12 we have obtained the following table:

\begin{matrix}n & v_{l,n}(2, 2) & v_{u,n}(2, 2) \\ 1 & 156.3722891740 & 156.5137544776 \\ 2 & 156.5105038907 & 156.5111569644 \\ 3 & 156.5111409414 & 156.5111441617 \\ 4 & 156.5111440813 & 156.5111440975 \\ 5 & 156.5111440970 & 156.5111440971 \end{matrix}

thus

156.5111440970 ≤ v_{l}(2, 2) ≤ v_{u}(2, 2) ≤ 156.5111440971.

References

[1] R. I. Avery, A closed form solution of a discrete correlated random walk, Journal of Difference Equations 12 (9) (2006), 949-956.

[2] R. I. Avery, A combinatorial approach to discrete diffusion on a segment, International Journal of Difference Equations 2 (2) (2007), 177-186.

[3] R. I. Avery, A note concerning discrete transportation, Advances in Dynamical Systems and Applications 3 (1) (2008), 25-28.

[4] R. I. Avery and G. Berman, A note concerning a discrete two dimensional diffusion problem and random walks, International Journal of Pure and Applied Mathematics 14 (1) (2004), 103-109.

[5] R. I. Avery and G. Berman, Discrete two dimensional diffusion problem on a half plane, Dynamic Systems and Applications 15 (2006), 35-42.

[6] C. de Boor, K. Höllig and S. Riemenschneider, Fundamental solutions for multivariate difference equations, Amer. J. Math. 111 (3) (1989), 403-415.

[7] S. S. Cheng, Partial Difference Equations, Taylor and Francis, New York, 2003.

[8] R. Courant, K. Friedrichs and H. Lewy, On partial difference equations of mathematical physics, IBM J. Res. Develop. 11 (1967), 215-234.

[9] D. Guo and V. Lakshmikantham, Nonlinear Problems in Abstract Cones, Academic Press, San Diego, 1988.

[10] W. G. Kelley and A. C. Peterson, Difference Equations: An Introduction with Applications, Academic Press, San Diego, 2001.

[11] K. H. Rosen, Handbook of Discrete and Combinatorial Mathematics, CRC Press, New York, 2000.

[12] J. Veit, Boundary value problems for partial difference equations, Multidimens. Systems Signal Process 7 (2) (1996), 113-134.

[13] B. Van Der Pol, Probability and Related Topics in Physical Sciences, Appendix IV: The Finite Difference Analogy of the Periodic Wave Equation and the Potential Equation, Interscience Publishers,Inc., New York, 1959.