Algorithms for the statistical modelling of the solution of elliptic-type boundary value problems

vo1.24,No.S,pp.153-157,1984

U.S.S.R. Comput.Maths.Math.Phys., Printed in Great Britain

CQ41-5553/84 $lO.OO+O.OO 01985 Pergamon Press Ltd.

ALGORITHMS FOR THE STATISTKAL PlODELLINGOF THE SOLUTION OF ELLIPTIC-TYPE BOUNDARY VALUE PROBLEMS* A.A. KRONBERG Methods of ptimizing random walk over spheres algorithms for the numeri.cal solution of the Dirichlet problem for the Laplace equation have been The finiteness of the variance of the estimate and the investigated. volume of computational work in the case of mixed boundary conditions are studied. Methods of statistical modelling which make use of random walks over spheres have been successfully applied to the solution of elliptic boundary value problems in special cases when the solution has to be found at a small number of points (or a small number of functionals or its are required to be estimated from the solution), and the region is multidimensional boundary is of fairly complex form /l/. Unlike many other areas of application of the various statistical modelling methods /2/, spheres algorithms have not been question relating to the optimization of random walk over In the present paper analogues of the principal-part isolation studied to any great extent. and essential-sampling methods /2/ for optimizing the numerical solution of the Dirichlet problem for the Laplace equation are investigated. In applications, boundary value problems are often encountered where a Dirichlet boundary condition is imposed on one part of the In this case boundary while a Neumann condition is imposed on the remainder of the boundary. algorithms associated with reflection in the boundary /3/ have been extensively used for a However, the finiteness of the dispersion in the estimates of the solution has long time. not been proved up to the present time and the question of the amount of computational work This omission will involved in a random walk over spheres scheme has not been investigated. be rectified below. Additionally, the problem for an inhomogeneous medium (containing dielectrics, for example), which has important practical applications at the present time, is considered. 1. The Dirichlet problem for the Laplace equation. Let D be a domain of an s-dimensional (s>2) Euclidean space with boudnary Dirichlet problem for the Laplace equation in the domain D is considered:

Au(z)=O, z=D,

r&=.1,

The

(1.1)

u(z)=(p(z), x=r, (P(Z)&.*.

To solve (l.l), we shall use a random walk over spheres algorithm. Let r.=(z~D:p(z, p(z, l?) is the distance from a point x to the boundary r of domain D. Let us r)
IMq.(zo)-44I
and is independent

of

E.

The solution

of (1.1) is representable

where n, is the external normal at the point y=r and We shall seek a solution of problem Dirichlet problem. u(x)=u(x)+

2

C&i

c(.z,y) is Green's (1.1) in the form

in the form:

function

for the

(xl ;

i-1

where the function constants We have

The variance

$(s)

are harmonic,

of the estimate

*Zh.vychisl.Mat.mat.Fiz.,

fairly simple, and previously

n.(.% C,,...,C.) for problem

24,10,1531-1537,1984 153

chosen and c, are unknown

(1.2) has the form

154

we shall select c, such that Dq.(r,) should be a minimum. lected we obtain a system of linear algebraic equations in c,: I

If the term

O(e)

is neg-

c)G(xc,Y)

UJ I_,

r

-----~,(Y)3,(Y)dY-lt,(z,)~j(l,)]c,= any

J--an,

aGc=o,Y)

(1.3)

g(Y)~,(Y)dY--cp(z,)rt;,(2,).

r

The fortunate choice of the functions significant reduction in the variance c, values.

*C(s) enables one, as shown in /4/, to obtain a very of the estimate of the solution, even with non-optimal

Proposition 1. If the functions q,(z) are linearly independent, differ from a constant, and no linear combination of the functions $8(Z) is a constant, then the determinant of the system of equations (1.3) differs from zero. Proof. Let us consider the linear envelope of the system of functions define the scalar product B of the functions (t,(z),cpl(z)es(@,(z)]:

The scalar product B defines in the case when q(y)+const

(see /5/, chapter

s an,

aG(xo,Y) -~2(y)dy-

r

(J

IX) a positive-definite

i%(z)1

Hermitian

and

metric

since,

~g4yldy)‘>O.

I)

i.e.

is a unitary space. The determinant is a Gram determinant which is constructed, {$,(x)) in the case of the functions $(z), from the scalar product B and the condition of proposition 1 is satisfied by virtue of the fact that a Gram determinant of linearly independent functions differ from zero /5/. The solution of the initial problem is divided up into two stages. In the first stage the matrix elements and the right-hand side of (1.3) are evaluated along the same trajectories. The system of equations (1.3) is then solved and, in the following stage, problem (1.2) is Since the terms cp(z,)$,(.r,) in (1.3) introduce an additional statistical error, in solved. certain cases, when this error is large, it is not possible to solve the minimization problem but the minimization problem M&'(zo, c,,...,c.),which leads /4/ to the system Dn.(zr, c,, . . ..c.) of equations

(1.4)

The determinant of system (1.4) differs from zero as a Gram determinant with respect to a conventional scalar product in the space L,(r). Let us now consider a method which is, in fact, an essential sampling method. Without loss of generality it may be assuemd that the solution U(s) of problem (1.1) is strictly positive and a harmonic function cp,(z)>&>O, which is strictly positive in Ii is known to us which behaves like U(Z) in the domain D. It is natural to organise the non-uniform modelling of the trajectories T on the surfaces of spheres using Neumann's method /2/ to simulate the next point according to the density c$~(s); z belongs to the surface of the corresponding sphere. Let us consider a random walk over spheres without absorption in r, (see /6/). In the case of uniform modelling the limiting points of the sequences (.$) belong to r with a probability of 1 and induce there a measure aG(~~,y)/an, (see /6/l. Let us assume that we solve the boundary value problem

A@,(s)=0 when s=D,

*o(z)=lpo(~) when

z~r

by modelling with respect to the density c%(s) * A probability measure may be introduced on the set of trajectories T for a random walk over spheres: a family of consistent measures is initially defined in finite intervals of a trajectory and a measure is subsequently extended to the whole set of trajectories using Kolmogorov's theorem. A measure generated by c@,,(z) in a set of sequences will be absolutely continuous with respect to the measure generated by uniform modelling and, consequently, the limiting points

of the sequence (bi)*. belong to r with a probability of 1 and, since the dispersion of the estimate qr. is equal to 0, the estimate of the solution of problem (1.5) itself is equal to $o(Zo) with a probability of 1. The followjng assertions follow from this. Proposition

2.

The limiting

points

aG(G,Y)

Proposition

3.

{g,),. induce

of the sequence

a measure

on r :

\Po(to)

The relationship

holds. The measure which is actually induced during the random walk over spheres will converge in the sense of a weak convergence of measures as E--O to the measure in Proposition 2. The question of the number of steps in the process during the random walk over c%(z) is interesting. If we return to the proof of the theorem in /7/, it is found that the proof follows in the present case with just a single change: in formula (2) from /7/ the number N will now be dependent, not only on the angle p, but also on the function qO(z). It is thereby shown that the asymptotic behaviour of the mean number of steps for a wide class of three dimensional domains is also not worse than Ilnel for a random walk using the essential sampling method. Various different combinations of the two methods described above are possible: by calculating c,, one may take

as (PO(x) and, inversely, the calculation is carried out in the first algorithm with a certain value of 9.(z). In /8/, the method of statistical modelling was applied as an auxiliary method to the solution of a boundary value problem by deterministic methods. The inverse situation, when qpo(s) is first found by a deterministic method with small demandsoncomputer time (see, for example, the methods in /9/b, is also possible.

2, Mixed

boundary

conditions.

1. Let a Dirichlet condition u(z)=q~(z) be prescribed on r,cT and a Neumann condition au(z)/&=w(z) be prescribed on r*==r\r, . Let us consider two versions of reflection from Tr: I) a reflection from the last point of a trajectory at a distance h>e in a direction opposite to the normal n,,;2) reflection from a point zr along the same direction also at a distance h. We assume that the corresponding intervals over which the reflection occurs lie entirely in the domain D. The following formulae hold:

uce~,=u(g,+h)+hw(E.)+O(h*), u(I,)-U(Z,+h)+hW(2~)+0(h*),

I u(zt)--u(EJ I
In both cases the estimate

A,=

y;) 21,

of the solution has the form (I.h('l)

12.lb) where zr is a point corresponding to the first ingress into the &-neighbourhood of r,, Z.,r is the number of reflecare the points on I'*where the i-th reflection occurs, and .&(T) tions. on account of the smallness of e the choice of this or that reflection method depends onthe,magnitudes of the first and second normal dertivatives of the solution on rl. The following proposition follows from (2.1). Proposition

4.

The relationship

ihW+O(e)

Mb(T)* h+O(h*)+O(e)

’

where Ap,(z) =O, tfD, pi(r)=O, z=r,. Sp,(z)ldn.=i, z=r, Whereupon it follows that Iu..~(~~)-u(zo)~~O(~)+O(~/~)+O(~). holds. Let us sum up. it is necessary to select eah;. the asymptotic To obtain an approximate SOlUtiOn,

amount

156

of computational work will be of the order of finiteness of the dispersion. Proposition

5.

Iln~[h-'. Let us study the question

of the

The relationship

Drl.,h(lO)~~'+2a)p,(~O)+2P_l+O(h?)+O(e),

holds. Lemma.

The inequality I

MS,,~(T)~~P-'~-'+O(~')TO(E) holds. The probability that, having been repelled from Proof. r,, the trajectory upon next entering r. will appear in the e-neighbourhood of rr, where the Dirichlet condition applies, p2(z.i-h)+O(e), depending on the choice of the method of rebound. is equal to p,(Er+h)+O(e) or We note that only a regular random walk over the surface of spheres is considered in Sect.2. For small values of h both of these probabilities are equal to 8P? (4 --h+O(hz)+O(e). dnzn

Let us replace

all the probabilities

considered

by

The majorant value is calculated in terms o,f a derivative function of the distribution and is proof of the lemma is obtained by a single coarsening equal to (2-hP)(hP)-',whereupon of The proof of Proposition 5 is obtained by the direct application of Proposition 4 and h. the lemma which has been proved to the expression Dq.,n(z~). It is characteristic that the majorant of the dispersion is independent of h. u(z)=~(z) be specified on the whole boundary r but 2. Let the Dirichlet condition the subdomain S, on the boundary of which the condition e,au/an--e&/an is specified, is located strictly within the domain D. The case, when there are several such subregions, is treated in a similar way to that described below. It is well known that, in this case of a random walk over spheres, trajectories which enter into the &-neighbourhood of the boundary S with probability el(e~+Ez)-' are repelled at a distance h inside subdomain S along the corresponding normal and, inside the subdomain D\S , with a probability er(e,i-e&-l. In this case the approximate solution is equal to MT(G), where z, is the point of the boundary r which is closest to the point where the trajectory enters r, for the first time. The variance of this estimate is obviously finite if ~~(z)IGlx+? Let us denote by M&,(T,S) the mean number of rebounds in the random walk over spheres trajectory with a start at a certain point belonging to the e-neighbourhood of the boundary of subdomain S. Let us introduce a function p,(s) into the treatment. This function &S is a solution of the boundary value problem Ap,(s)=o,~~D\S,p~(z)=o on the boundary S, p3(z)= Let us put I, te.

It is then easy to obtain

an upper estimate

for the mean number

M%.,r(T, S)G(e,+ez)[e,(P,h+O(h*)+O(E))]-‘. then the mean number of steps will be asymptotically not worse than it is not worse than [ps(z,)f(l-pr(~,))XIM~.,h(T, S)] IIn El. s)+l] jln el; if x,=D\S, If (this case will be considered below) of the points of reflection is regular on as formula If

Z&s,

[M%.,h(T, the

measure

then the

(2.2) holds. As in the case of paragraph 1, it is necessary to choose eO, it makes sense to talk of a stationary distribution

157

Then the random walk over spheres algorithm, in terms of psi using the well-known method. e-neighbourhood of &S from within S, will be "repelled" into a random on entering into the point of the subdomain D\S which lies at a distance h from the central point of the subdomain S, where k is a random number with a distribution density y*. Fortunately, in applica(s=2) or a sphere (s>2). In this tions, the subdomain S may often be assumed to be circular casenopreliminarywork is required and a point is modelled on surface S with a regular density from which repulsion occurs at a distance h into the subdomain D\S. On account of this, when s=2. it is possible to pass to the solution of an analogous problem when the form of s is a unit sphere by means of the well-known conformal mapping of s into a unit by virtue of the invariance of the Laplace operator with respect to conformal circle IzIG mappings and the conservation of angles. Since infinitely many such conformal mappings exist, the question arises of theiroptimal choice (in this case optimality has the meaning of the smallest number of steps in the By virtue of equation (2.2), it is necessary to maximize the value of the trajectories). intearal

The solution of the optimization problem is elementary: invariant to conformal transformations /lo, page 80/.

the value of the integral

is

REFERENCES ELEPOV B.S. et al., The solution of boundary value problems by the Monte-Carlo method (Reshenie kraevykh zadach metodom Monte-Karlo), Nauka, Novosibirsk, 1980. and MIKBAILOV G-A., A course in statistical modelling (Kurs statisticheskogo 2. ERMAKOV S.M. modelirovaniya), Nauka, Moscow, 1376. (tletod statisticheskikh 3. BUSLENKO N.P. et al., The method of statistical experiments ispyanii), Fizmatgiz, Moscow, 1962. 4. KRONBERG A.A., On several algorithms for solving elliptic boundary value problems, Zh. vychisl. Mat. mat. Fiz., 22, 6, 1401-1409, 1982. 1966. 5. GANTMAKRER F.R., The theory of matrices (Teoriya matrits), Nauka, Moscow, Some continuous Monte-Carlo methods for the Dirichlet problem, Ann. Math. 6. MULLER M.E., Statistics. 27, 3, 569-589, 1956. e7. KRONBERG A.A., On the asymptotic form of the expectation of the number of steps in an spherical process (Ob asimptotike matematicheskogo ozhidaniya chisla shagov e-sfericheskogo protsessa), Zh. vychisl. Mat. mat. Fiz., 20, 2, 528-531, 1980. 8. ELEPOV B.S. and IL'IN V-P., On the combination of difference algorithms with the Monte-Carlo method, Zh. vychisl. Mat. mat. Fizz., 15, 2, 384-390, 1975. 9. ALEKSIDZE M-A., The solution of boundary value problems by expansion in terms of nonorthogonal functions (Reshenie granichnykh zadach metodom razlozheniya po neortogonal'nym funktsiyam), Nauka, Moscow, 1978. .. 10. POLYA G. and SZEGO G., Isoperimetric inequalities in mathematical physics, Russian translation, Fizmatgiz, Moscow, 1962. 1.

Translated

by E.L.S.