The efficiency of the use of asymptotic solutions in calculations by the Monte Carlo method

THE EFFICIENCY OF THE USE OF ASYMPTOTIC SOLUTIONS IN CALCULATIONS BY THE MONTE CARLO METHOD* B. A. KARGIN and G. A. MIKHAILOV Novosibirsk

(Received

7 September

1970)

VARIOUS algorithms of the Monte Carlo method, constructed solutions

of the adjoint problem, are known [l-3].

by using approximate

In the present paper the effi-

ciency of some algorithms of this type is studied by the example of the passage of particles

through a plane layer of matter with isotropic

scattering.

Algorithms

constructed

by means of the approximate asymptotic solution of Milne’s spherical

problem for anisotropic scattering are also considered. These algorithms are used to solve practical problems of estimating the time distribution of the illumination of an area at a great optical distance from a directed source of light.

1. Introduction The process of particle

transfer can be described

by an integral equation

of the second kind [41

where f(x) is the density of collisions

(or of scattered

particles)

in the phase

are usually calculated space X. Linear functionals of the form 1, = (f, q) by the Monte Carlo method. The variance of the estimates can be decreased [l, 21 by using information about the solution of the adjoint equation f’ =

K’f’

+ cpLet

g(z)

> 0

probability density g(z)

/ (q,, g)

and (5,) k (z’, z) g(z)

be a homogeneous Markov chain with transition / [K’g]

(z’),

and absorption probability n,(x).

Qn by the relations *Z/z. uychisl. Mat. mat. Fiz.,

12, 1. 150-158, 1972.

187

initial density \c,(5) X We define the random weight

B. A. Kargin and G. A. Mikhailou

188

vo=- (44g)

=

Qn

g(xo) ’

It is known [51, that is p,,(z)

>

Qn-i

rK'g1

h-i)

g(xn) [I -

Pn-1

(Gl-1)] *

1, N

(2)

I. = ME,

me

z

EI

Qndxn).

n=o

Here N is the ordinal that if p(x) >

(3)

number of the last state

It is also known

0 for 4(x) > 0, rEe ,,=

J.=h,

Three algorithms described

of the chain.

Q*b) P&v) -

of the Monte Carlo method for estimating

the functional

J9 are

below.

Algorithm

A.

This

algorithm

is based

on relation

(2) with the condition

p, I 0, if n& m, and p, > 6 > 0 if n > m. It is known (see [ll), that the variance DE for this algorithm equals zero for g = fr and m = 00. Some estimates of the variance

tor the algorithm

Algorithm Therefore,

B.

This

algorithm

the additional

in [2l that

A when g = fr were obtained

DE < 1;2

is based

in [ll.

on the relation

(2) for p =

condition Kg/g < i is required for this algorithm if g = f”.

here.

1-

Algorithm C. This algorithm is based on relation (3) for p = 1 If g = p the variance of this algorithm Dn = 0 (see, for example, 121). In the next section we give the results of calculations particles through a plane layer of matter. The asymptotic problem

is used es the function

of the algorithms estimate

is also easily

obtained

computer

time required

to attain

the transmission

g(x).

A and B for this case

probability,

Theoretical were

for the Algorithm a specified

increases

relative

in proportion

C.

K*g/g

3

on the passage of solution of Milne’s

estimates

obtained

K*g/g.

It was shown

of the efficiency

in [ll and 131. A similar It is found that the accuracy

in estimating

to the optical

thickness

the layer, and not exponentially, as in the direct simulation of the partical transfer process. The theoretical results obtained are not precise and do not show that any of the algorithms considered has any advantage over the others. Hence a numerical comparison of them with one another and with direct simulation is of interest. However, we notice that the algorithm A has an obvious advantage: its realization does not require the condition K*g/g f 1, which

of

0.

~fficje~cy

is not satisfied,

of the use of asymptotic

189

solutions

for example, if g is an approximate asymptotic solution.

2. Numerical

comparison

of the various

algorithms

We consider ah omogeneous layer of matter 0 < z < 20 in which the mean free path of the particles is E = 1, the absorption probability at a collision

P =I-q=

0.2 and the scattering

is isotropic. A beam of particles falls on the surface z = 0 in the direction of the z-axis; the angular distribution of the particles emitted through the surface z = 20 has to be estimated by construe ting a histogram with step Ap = 0.1 of the cosine of the angle between the direc tion of motion of a particle and the axis. Therefore, the problem consists of estimating the probability J,

of the passage

of particles

with the latter direction

in the interval

O.l(k-- 1) \( p < Let f(z, IL) be the particle

k = 1,. . . , 10.

O.lk,

scattering

density.

It is known 161, that

J,+= (f, tpk),

where

The asymptotic solution (that is, the solution in the depth of an optically

thick

layer) of the adjoint equation for this problem, is independent of the number k and is proportional to the function [?‘I:

g(z, j_k)= erlL(l -

(4)

@/I;) -’ = e*lLu ( p) .

The diffusion length L is defined by the equation

1 /q =

(L / 21) In [(L +

I) /(L-Z)]. It was shown in Ill that the simulation of the Markov chain corresponding to the expression (4) consists of the following: the length of the mean free path of the particles

is simulated as for the “exponential

with parameter l/L “without emergence”

transformation”

[81, and the direction of scattering

is simulated according to the density proportional to a(~>. Using the explicit form of the operator K*, it is easy to obtain ill, that here

Wgl(z, PI’---&?(z, PI-- a(cl)b(z*CL),

(5) where

b (z, p) =exp(--z/Zlpl)

Expression

if&CO, and exp{-(z,-z)/Zp} ifp>O. (5) shows that K*g/g & 1, and therefore the algorithms B and C

can be used in this case.

B. A. Kargin and

190

G.A. Mikhailov

TABLE 1 Alm=50.

I

-

i-

7)

t=223)

wt

-c ft=57.4

f .$ . lO-*mk.iO" -t

0.8 0.9 144 199 1.0 ,282 Table trajectories

3.21 17.6 250 11.3 7.36 43.7 23.1 6.51 34.2 14.8 40.5 9.85 55.6 ";z '78.8 2:67 5.15 111 2.17 3.80 150 i.56 1.96 200 1.39 1.56 286 1.16 1.09

b% 0:52

1 shows the results of the calculations. were simulated. The following notation

tion of the calculation

(in seconds)

I

-I

2 tajc40-

"1( %

4.45 15.9 145 ii.6 14.1 114 9.58 52.7 $2.4 36.7 5.83 19.5 5.36 16.5 52.6 3.57 79.3 7.29 3.76 !% :2; ;-z 203 1:41 0.95 1:29 262

For each algorithm 10 000 is adopted: t is the dura-

on the BESM-6 computer,

mk is the estimate

of the probability Jk, and ok is the corresponding estimate of the relative error probability (as a percentage). The quantity to* is a criterion for estimating the quality

of the algorithm.

For the algorithm

A the quantity

by starting with the fact that the sum of the series ciently small in comparison with the quantity

The result

of the preliminary

estimate

m was first estimated

OD (K”% c n=m

was as follows:

VA)

m = 60.

is suffi-

Calculations

showed that the optimal value of m lies between 50 and 60, and that the efficiency of the algorithm changes little when the value of m is varied in this interval. The value

m = 50 was used in the calculations.

absorption

was simulated

Analysis

with probability

of the results

showed

For numbers

of collisions

p = 0.2 in the algorithm

that here the algorithms

n > SO the

A.

A and B have

practically the same efficiency and are rather better in quality than the algorithm C. In general the experiment may be considered to have ended in favour of the algorithm rigorous

A, since condition

the realization K’g/g

4

of B and C requires

the satisfaction

of the

1.

Additional analysis showed that the computer time expended here was approximately only 10 -’ of the time required to attain the same accuracy by direct simulation of the particle transfer process.

Efficiency

of the use of asymptotic solutions

3. Connection with the exponential It is well known on the representation

[81that

191

transformation

the algorithms of the Monte Carlo method are based

of the beam of particles

@(z, cl> in the form

where c < u is some constant, u = l/Z. The integrodifferential equation for @ii obtained in this way may be regarded as the equation of particle transfer in a fictitious

medium for which (O - cp) is the total cross-section

and os/(u - C,Y)

is the survival probability; us = uq. A transformation of this kind for spherical geometry (that is, when the transmission of particles to a point detector is considered) was studied in [91. If the quantity qi = (~~/(a - cp) > 1, it is necessary

either to simulate the multiplication,

or to multiply the weight Q by ql,

which may lead to a considerable increase in the variance. An attempt may be made to compensate for this fictitious ‘weight” by a special modification of the scattering: where a(p) is an as yet unknown @JIL’, cl)= 4l4 P’)4:l& function. After the simulation (without absorption) of such a scattering the “weight” is transformed by the formula

It is obvious from this that if the relation +? (1-

cZp’)a(p’)

= q J a(~, v’)a(+p, --i

is satisfied, the fluctuations of the weights do not accumulate. Thus we have arrived at the characteristic equation of Milne’s problem [71 by improving the “exponential” transformation, without using the asymptotic solution of the adjoint equation.

4. Use of the approximate asymptotic solution of Milne’s spherical problem for anisotropic scattering Let the layer

0 < z <

homogeneous substance, following notation:

r =

H of three-dimensional space be filled with a absorbing and scattering particles. We introduce the (z -

5*,

y-

y*,

z-

z’)

is the radius vector of the

point (x, y, z) relative to the point (z*, y’, z’), 1r 1 = r, co’ is the direction of motion of a particle before a collision, ois the direction of motion of a particle after a collision, p = UO’ is the cosine of the scattering angle. The optical properties of the medium are defined by the extinction coefficient (a

192

B. A. Kargin

completely

efficient

cross

and G. A. Mikhailov

section)

2, the survival

indicatrix

w(p), which is the probability

scattering

angle.

We assume required

that the particle

to determine

source

the distribution

density

q

probability distribution

is situated

of the cosine

on the plane

of the particle

and the scattering

x = 0;

it is

flux in some region

of

space localized about the point (x*, y*, z*), which is at a great optical from the source. For the solution of Milne’s

adjoint

of this problem spherical

it is natural

problem

[7].

distance

to use the asymptotic

Apart from a constant

factor,

of the

solution this solu-

tion has the form

g(r, pL,)= exp(--r/~)a(k)lr, where

p, =

connected

and the ‘diffusion

o(r/r),

length”

L and the function

a(/.+) are

by the equation

(6) A simple

investigation

shows

that because

tor the use of g(r, p,) is extremely

g0(r,

integral

transport

of the Markov chain

equation

corresponding

The mean free path h is simulated tional

Hence

of r in the denomina-

the approximate

function

form of the kernel

of the

cl?)= exp(--ri-jL)a(p,). Using

was used to modify the simulation. general

of the presence

inconvenient.

the explicit

[41. it is easy to obtain to expression according

(7) consists

to the density,

that the simulation of the following. which is propor-

to the function

where

(8) Here

r(Z) = 7/[P(O)+ Z2+ 2WOb)l. r(Z)=(z-z*,

y--g’,

and r(l)=r(O)+ lo, r(Z)=(r(Z) we have

Z-Z’), I. Since

(x, y, z) is the point of the next collision, p,(Z)=or(Z)/r(Z)

=

ar(Z)/dZ,

Efficiency

of the use of asymptotic

193

solutions

F(l) = j kb(t)dt + r(O).

(9)

0 Consequently,

This means particle is,the

apart from a multiplier

that the Markov chain

transport

process

total extinction

the length

corresponding

in a fictitious

coefficient.

It follows

where Q is a random number uniformly of equation

(11) can easily

displacement Q(l)=exp

The same algorithm for the “spherical a parameter

l/L

a

(Z -I- ~1, /L) the equation

distributed

in the segment by using

by the Gweighting”

IO, 11. The

the expression

(8).

The

factor:

(““‘~““‘)/(i+Jf$-). for simulating

completely

The scattering

to solve

be obtained

is compensated

exponential

tion of the scattering

represents

pT),

Ina

l+(F(l)-F(O))/~~=-P-,

resulting

go(F,

from (9) and (10) that to simulate

1 of the mean free path it is necessary

(11) solution

to

medium for which

the length

of the mean free path is obtained

transformation” determined

is carried

I9I.

by equation

out in agreement

must be simulated

according

The difference (6) is used,

is that here and the simula-

with this modification. to the density,

proportional

to

the function

(12)

h(a)=

wb)4P*).

Here the Uweight” of a particle

must be multiplied

by the quantity

where

The algorithm is considerably simplified if the asymptotic solution is found in the transport approximation [lOI, which leads to isotropic scattering with the

t?. A. Karginand G. A. Mikhailov

194

effective

mean free path

and effective

survival probability - Q-

Pi1- PI i--@m

,.where

p =

* I+lpW(lr)d~ --i

Then L is the solution of the transcendental

equation [71

and -1

a(pr)=(l+g) In this case

5. The method described

A numerical

example

was used to solve the problem of estimating

the time

distribution of the illumination of a small area situated at a large optical distance

from a directed source of light.

uniformly over the area of the circle tropically

in a cone of directions

The particles

by the density, proportional to the function u = 0.023.

time).

ye+,

coefficient

determine the time distribution

The

from the source is defined where

The velocity of propagation of particles

The extinction

and iso-

with aperture angle 30” about the x-axis.

distribution with respect to the time y of the particles dimensionless

(photons) are emitted

t = 0, y* + Z* Q 25 . fO-*

1~= 9 (y is the in the medium

2 = 1, and q = 0.9. It is required to

of the beam of particles

arriving at a circular

area of radius 0.1 with centre at the point (H, 0, 0), oriented perpendicular to the x-axis (that is, the time distribution particles

intersecting

we took the scattering

of illumination of the area).

the area are counted as absorbed. indicatrix

of sea water (II= 0.9).

table with linear interpolation between the nodes.

The

For the calculations It was given as a

In this case the integral for

Efficiency

TABLE

of the use of asymptotic

TABLE

2

195

solutions

3

Ek.lO’

-i-l-l 24.8 0 1.0 i 0.602

I$‘)

is expressed

&+#I

+

0.296 0.947 1.260 0.984 0.831 0.612 0.440 0.242 0.222 0.210 0.189 0.108

20.70 21.15 21.60 22.05 22.50 22.95 23.40 23.85 24.30 24.75 25.20 25.65

23.0 1.003 0.578

in terms of elementary

120

2.82

:ti 119 183 201 227 289 306 365 446 693

00~~ 0:85 x I:89 3.69 3.71 4.00 4.41 7.71

functions

pn+rlZL+ bi + ~17’

n+r)ZL1n ( p,IlZL +ba +p,’

pr’( w, - w

where pr’ # -4, pr’ # 1, b

b.=~l~~r.+r~')'+[(zL)~-*](i-p,")), N is the number of nodes

at which the indicatrix

is specified;

N-l

w.b+i - 1

(WUn+,Pn -

ZLpn -

1

(l-h- Pn+i)-I,

N-l

Z(1)=

The simulation

w -

1)

of p according

c [(WI-

Wn+i)(Pn+i - pn)ZL +

to the density

(12) can be carried

out by the method

B. A. Kurgin

196

of elimination according bility

[ll],

choosing

to the function

of elimination

p according

to w(u) and carrying

(1 + l,t, / ZL)-‘.

is large,

with the approximate

and G. A. Mikhailov

the simulation

For values

out the elimination

of pr where the proba-

of ~1can be carried

out in accordance

density:

where

As a check ancy”

on the applicability

of equation

of the transport

(6) was calculated.

Table

approximation

the “discrep-

2 gives the values of the left and It is seen from Table 2 that the two

right sides of equation (6) for three points. sides of the equation differ by not more than 7%. The calculations 20, 30.

Table

illumination described

were carried

3 shows the results

of calculations

distances

with direct

simulation.

the time nodes of the histogram oPk

The following

notation

(in dimensionless

and ask are the relative

of E, (in percentages), tP = 300 respectively;

‘G=

units),

It is seen from Table computer

was adopted:

of the

yk are

E, is an estimate

probabilities

of

of error in the

for direct simulation and using the asymptotic and t, = 1100 are the times for the calcula-

3 that in order to attain described

time that is required

10,

of the algorithm

solution, tions of E (in seconds) on the BESM-6 computer (36 000 trajectories), simulation and using the asymptotic solution, respectively.

tion, the use of the algorithm

XH =

of the time distribution

E of the area for I = 20, and also a comparison

the illumination, calculation

out for three optical

for direct

requires

a given accuracy about lo-’ times

for direct

of calculathe amount of

simulation. Translated

by J. Berry

REFERENCES 1.

MIKHAILOV,

G. A.

The use of approximate

improve the algorithm

solutions

of the conjugate

of the Monte Carlo method. Zh. uj;chisl.

problem to

Mat. mat. Fiz.,

9, 5, 1145-1152, 1969. 2.

COVEYON, R. R., CAIN, V. R. and YOST, K. J. Adjoint and impotance Carlo application. Nucl. Sci. Engng, 27, 2, 219-234, 1967.

in Monte

3.

MIKHAILOV, G. A. The use of asymptotic solutions of Milne’s problem in calculations by the Monte Carlo method. Izu. VUZou. Ser. fiz., 10, 76-79, 1968.

Efficiency

4.

ERMAKOV,

of the use of asymptotic

S. M. and ZOLOTYKHIN,

to calculate

Shielding

protection

(Vopr.

Reactors

V. G.

Application

from nuclear radiation, fiz.

zashchity

197

solutions

of the Monte Carlo method

in: Problems

reaktorov),

of the Physics

149-174, Atomizdat,

of Moscow

1963.

5.

I
D.

Monte Carlo methods for the iteration

of linear operators.

Usp. mat.

Nauk, B2, 5, 149-174. 1957. 6.

BUSLENKO,

N. P. et al.

The Method of Statistical

(Meted statisticheskikh

7.

DAVISON,

B.

Moscow,

8.

FANO,

Transport

U.

Gamma Radiation

MIKHAILOV,

G. A.

Transport

Application

Method in the Problem probleme perenos MARCHUK,

G. I.

nykh reaktorov). 11.

BUTLER,

Theory (Teoriya

Fizmatgiz,

perenos neitronov).

Moscow,

1962.

Atomizdat,

(Perenos

gamma-izlucheniya.)

Atomizdat,

1963.

of the transport of particles

10.

(Monte Carlo method)

Trials

(meted Monte-Karlo)).

1960.

Moscow, 9.

Neutron

ispytanii

of the Transport

izluchenii),

Methods

Moscow,

sampling

sium on Monte Carlo Methods,

transformation

ofRadiation

67-71, Atomizdat,

of Designing

Atomizdat,

J. W. Machine

of the exponential

by the Monte Carlo method.

in calculations

In: The Monte Carlo

(Meted Monte-Karl0

Moscow,

Nuclear Reactors

v

1967. (Metody rascheta

yader-

1961.

from given probability

distributions,

249-264, Wiley, New York, 1956.

in: Sympo-