by DR. EVERETT L. JOHNSON The Boeing Company
Wichita, Kansas
INTRODUCTION
The work done to date on the analog/hybrid Monte Carlo solutions of partial differential equations can be summarized by reviewing three works: one research report and two Ph.D. theses.
Chuang, Kazda, and Windenecht were the first to demonstrate the feasibility of a Monte Carlo solution of a class of partial differential equations on an analog computer.! The boundary value problems for which their stochastic solution technique is applicable belong to a family of generalized Dirichlet problems of the form
(1) where Kl (Xl, X2) and K2 (Xl, X2) are arbitrary functions of Xl and X2. The boundary, c, is an arbitrary, finite closed curve-a Jordan curve.
The boundary-value function ¢(XI, X2) is a bounded, single-valued piecewise continuous function of Xl and
X2. DI and D2 are constants.
The method developed is based on the direct relation that exists between partial differential equations and the random process that arises in the analysis of electric circuits subjected to random excitations.2
The electrical equations to be simulated for the solution of Equation 1 are
dXI
at +
KI(XI, X2) = NI(t)(2)
where NI(t) and N2(t) are noise generators with Gaus-19
sian amplitude distribution and power spectral densities
D! and D2 respectively. Figure 1 demonstrates the simulation technique.
Fundamental to the research of Chuang, et aI, is the theorem of Petrowsky3 which guarantees the con-vergence of the stochastically obtained solution to that of the generalized Dirichlet problem, Equation 1.
The solutions to several one and two dimensional problems were given in the paper by Chuang, et al.
The solutions were obtained as follows: An analog program as shown in Figure 1 was patched. The initial conditions of integrators 1 and 2 were set to the co-ordinates of the point for which the solution was de-sired. Boundary crossings were detected by using an oscilloscope, bounded region mask, and photo tube.
Upon detecting a boundary crossing, the value of the boundary-value function at the point of intersection was recorded and the process repeated. The approxi-mate solution was then given by
1 N
f(~) = N
E
¢(Sn) (3)where Sn is the coordinate of the nth crossing, ~ the solution point, and N is the number of repetitions of the procedure. The simulation equipment allowed up to 2000 random walks per hour. Errors were in the .5 percent range and were attributed primarily to sta-tistical variations, presumably in the noise source.
Little4 extended the class of partial differential equa-tions for which the methods developed by Chuang, et aI, are applicable while developing a technique of solution utilizing an analog computer linked toa digital computer. Little solved three types of partial differ-ential equations: parabolic, elliptic, and non-homo-geneous.
Little used the analog computer for simulating an electric circuit excited by random noise. The digital
20 Spring Joint Computer Conference, 1970
BOUNDARY DETECTION
ANALOG - - - - MODE
CONTROL LOGIC
~=gi I--_~ AVERAGING
GENERATOR DEVICE
Figure I-Continuous random walk program
computer collected and averaged the resulting bound-ary values and controlled the modes of the analog components. With his hybrid system, EAI 231R-V analog computer and Logistics Research Alwac III-E digital computer, he was able to obtain 10 random walks per second. The Monte Carlo solutions compared favorably with the analytical solutions of the example problems.
Handler5 demonstrated that by use of a high speed repetitive operation analog computer with parallel logic capability, the ASTRAC II, that the Monte Carlo solution techniques developed by Chuang, et aI, and Little could be made competitive with more conven-tional finite difference digital solution techniques. In fact, he demonstrated the ability to plot continuously and directly the solution to partial differential equa-tions. This was accoIllplished by slowly changing the coordinates of the point for which the solution was desired while performing 1000 random walks per second.
The averaging of the intersected boundary values was done by a simple analog averaging circuit.
A. W. J\!Iarshall6 has stated that if a random sampling method is to be used to solve' a problem, attention should be turned to three topics.
(1) Choosing or modeling the probability process to be sampled (in some cases this means choice of the analog; in others a choice between alternative prob-ability models of the same process).
(2) Deciding how to generate random variables from some given probability distributions in an efficient way.
(3) Variance reduction techniques, i.e., ways of in-creasing the efficiency of the estimates obtained from the sampling process.
In summary:
(1) Petrowsky has defined a class of stochastic proc-esses which can be used for obtaining a solution to the Dirichlet problem.
(2) Wang, et aI, have established a random process, in the class defined by Petrowsky, to be sampled.
The first two topics suggested by Marshall have been considered previously and satisfactory results obtained.
The third topic is the subject of this paper.
The means toward the end will be an examination and implementation of a technique of stratified sam-pling using the Green's function for a rectangle and for a circle.
The technique is implemented by choosing an ap-proximating region, Ra , which is totally contained within the region, R, for which a solution is desired. Portions of the boundary of Ra may be coincident with the boundary of R. For a point contained in Ra the solution is found by performing walks which originate from the boundary of Ra. The number of walks, N m, which originate from the mth segment is given by the negative of the normal derivative of the Green's function for Ra integrated over the mth segment. A substantial re-duction in variance results from the use of the technique.
A VARIANCE REDUCTION TECHNIQUE FOR .THE CONTINUOUS RANDOM WALK For a solution of Laplace's equation for the region R of Figure 2, the continuous walk technique moves
I
,
I, ,
-,
-\
I
\
"
Figure 2-Region of solution with approximating region
continuously from the point of interest, ~, in the region R until a boundary is intersected. The value of the boundary function cP (s) is recorded and after N walks the solution is
(4) If the Green's function for Laplace's equation and the region R are known then the solution can be written7 as
f
aGcp(~) = - cp(s) - (~, s) ds
c an (5)
where ~ represents the coordinates of the point for which a solution is desired, s is the variable on c the boundary of the region R, and (aGjan) (~, s) is the derivative of the Green's function with respect to the normal vector of the boundary of R.
In Figure 2 consider the region Ra contained by R with the boundary represented by a dotted line. Each walk leaving the point ~ must intersect the boundary of Ra at least once before intersecting the boundary of R. If a boundary function CPa (Sa) were given for Ra, where Sa is the variable on the boundary of Ra, then CPa (Sa) evaluated at the points of first intersection with the boundary of Ra allows the solution for Laplace's equation in Ra to be written
1 N
CPa(~) =
N E
cpa (Sai) . (6) It was reasoned that if Ra were a region with' a known solution this information might be used to determine the points· of intersection on the boundary of Ra for walks originating at ~. Due to the Markovian nature of the random walk process, walks originating from the boundary of Ra with the correct distribution for con-tinuous walks from ~ should give the solution to Laplace's equation for R. The information sought can be obtained from the normal derivative of the Green's function for Laplace's equation in Ra. The probability density function properties of the Green's function are well established.8 A heuristic argument is given in Appendix A to support the use of the Green's function to find the proportion of the N walks with origin at ~which intersect any segment of the boundary of Ra.
If the boundary intersection coordinates are stored during a random walk solution of Laplace's equation for the region R then these values can be used for the solution to Laplace's equation for various boundary functions. The set of coordinate points is then' an ap-proximation to the Green's function for the region R and the walk origin ~. It is important that these inter-section coordinates have as nearly as possible the same statistical parameters as the Green's function. The
im-Variance Reduction Technique 21
portance is obvious if the boundary function IS ex-pressed as a power series.
The solution to Laplace's equation is then seen to be a function of the moments of S
cp(~)
=f
(ao+
alS+
a2s2+ ...
)Gn(s) ds. (8)8
The possibility of obtaining a variance reduction in the sample average by originating continuous walks from theoretically determined origins on the boundary of Ra prompted the research reported in this paper.
A STRATIFIED SAMPLING TECHNIQUE USING GREEN'S FUNCTIONS
Stratified sampling is a sampling technique which gives a reduction in sample variance if the population from which samples are to be made can be divided into sub groups which have variances smaller than the original population. The sub group selection and the number of samples to be made from each sub group must be selected such that the parameters to be de-termined from the samples are the same for the stratified sampling as for samples taken from the original popu-lation. The use of Green's functions to determin.e the distribution of walk origins on an approximating region Ra contained by R permits the division of the original sample popUlation that exists at ~ to be divided into sub groups with variances smaller than that of the population at ~. The technique described below was implemented and shown to give a variance reduction before its recognition as an application of stratified sampling.
The technique consists of performing walks from the boundary of the approximating region Ra which con-tains the point, ~, for which a solution is sought. The boundary of Ra is divided into M segments, ~Ai' The proportion, pm, of the total number of walks, N, to be originated from the mth segment of Ra is determined by (9) where Gna(~, Sa) is the normal derivative of the Green's function for Ra. Note that pm is the probability of a walk originating at ~ intersecting the boundary segment
~Am' The segments must be small enough that the solution to Laplace's equation at the mid points of two adjacent segments does not differ greatly.
Consider the solution of Laplace's equation using the technique described above for a region R, Green's
22 Spring Joint Computer Conference, 1970
function with normal derivative Gn(~, s) and with boundary function c/>(s) = s on a portion of the bound-ary, ~s, and zero elsewhere. Divide the boundary of the approximating region into M equal segments ~Am.
For the problem described, the solution is 1 N
S = - I:Sk N k=l
(10) where
S
is the average value of the intersection co-ordinates on ~S. The variance of the solution is9 with Ra1 M
YeS) = -
I:
PmO"m2N m=l (11)
without Ra
giving a variance reduction
o(V) =
~(
Qu2+
2f:
Pm(Um - U)2 -f:
Qmu",'Pm) (13) wherepm = the probability of intersecting ~Am
O"m2 = the variance of the average value of inter-section coordinates from walks originating at
~m, the mid point of ~Am
Um = the expected value of the intersection coordi-nates resulting from walks originating at ~m U = the expected value of the intersection
coordi-nates resulting from walks originating' at ~
Q = -
fGn(~'
s) ds8
Qm = -
f Gn(~m,
s) ds.8
The wavy line under
S
in Equation 12 denotes the fact that in order to get the equation in terms of Equation 11 it was necessary to assume that any walk originating from ~ and intersecting a segment of the boundary of Ra continued its motion from the mid point of the segment. For sufficiently large M the difference in results is negligible. For the experimental verification given below a value of M = N was used.For verification of the above results, Laplace's
equa-7 ~(l.X) - 0
1 r---~~~----~----~
o~---~-(-o-.x-)-.-. -X---~l x
Figure 3
tion was solved for the region and boundary values given in Figure 3.
A circle with a radius of .45 centered at (.5, .5) was used for Ra. Solutions were obtained for x = y = .5.
Approximately 100 walks per second were made-this excludes the time required to compute the number of walks made from each segment.
For each value of N, ten solutions were made. The ten values were used to compute the standard deviation of the sample average. The sample average is the
so-2S
C"'\
.... o H o
20
15
10
5
\
\
, ,
, , , ,
, , ,
--- Theoret1oal W1thout Ba
_ _ Theoret1oal W1 th Ra
' ... , ...
...
--- --- ---
---+
o L---4TO-0--~-8~0"0---lT20-0----~1~6~OO~--~200'0
.UIIl~r or Walke
Figure 4-Standard deviation for X = .5, Y = ~5, with Ra
lution to the given problem. The experimental values were plotted along with theoretical values for (J' for the cases with and without Ra using the square root of Equations 11 and 12.
The problem was programmed on an EAr 690 Hybrid Computer. The digital portion of the hybrid system computed the initial coordinate values for the walks, the number of walks to be made from a segment, re-corded the boundary intersection coordinates, com-puted the average, and controlled the modes of the analog computer. The analog portion integrated the noise to generate the walk, detected boundary improve-ment in standard deviation using Ra is apparent.
By using Ra for the solution of Laplace's equation the following benefit is realized:
for x = y = .5, a ,158% reduction in variance or the same variance as 'without Ra with 42% as many walks.
For the problem of Figure 3, Figure 6 gives the
theo-60 that for a given variance in sample average the number of \valks required is not a constant.
The technique described and demonstrated in this section gives a substantial reduction in error in the continuous random walk solution of Laplace's equation.
For a given number of walks the amount of error reduction is dependent on the position of the point for which a solution is sought. The next section contains an extension to the technique which in many cases Figure 7-Variance with Ra as a function of radius of circular Ra
centered at. X - Y = .;") for t.he point X = Y = .• 5
24 Spring Joint Computer Conference, 1970
VARIANCE IMPROVEMENT USING BOUNDARY COINCIDENCE
The continuous random walk solution of Laplace's equation with boundary function q, (s) for a region R with approximating region Ra can be written as
M Nm 1 N m
q,(~) =
EN
NmE
q,(Sim) (14)or
M
q,(~) =
L
Pmq,(~m). (15)m=l
q,(~m) is the sample average at the mid point of the mth segment, N m is the number of walks intersecting the mth segment of Ra and Sim is the coordinate of the intersection of the ith walk from the mth segment.
Equation 15 is the continuous random walk solution of Laplace's equation in Ra with boundary function
q,(~m) on ~Am. Since pm may be computed using the approximating region Green's function, Equation 15 can be written
(16) or
= -
i2 f </>(~m)Gna(~,
s) dsm=l ~A,,~
(17) where Gna is the normal derivative of the Green's function for Ra. If Ra can be positioned such that portions of its boundary are coincident with the bound-ary of R then the boundary function on the coincident portions of Ra need not be approximated by q,(~m) but is the same as the boundary function for R. Equation 17 can then be written
q,(~)
=-! q,(S)Gna(~,
s) ds+ t Pmq,(~m)
(18)8e m=l
where Se is the coincident portion of the boundaries and L is the number of segments not coincident.
The second term of Equation 18 can be written
which can be written 1 K
- L:
q,(Si)N i=l
where Si is the ith intersection coordinate,
L
K =
LN
m ,m=l
(19)
(20)
(21)
and
0< K ~ N. (22)
K is the number of walks actually made and is de-pendent on the extent of the boundary coincidence and the coordinates of the point for which a solution is sought. N will now be referred to as the base number of walks.
Equation 11 gives the variance when Ra is used which is the weighted sum of the variance from the M segments of Ra. The variance of the coincident seg-ments is zero giving an additional decrease in the variance.
EXAMPLE PROBLEM
The solution to Laplace's equation is desired for the region shown in Figure 8 and boundary conditions:
q,(0, B) = (3
</>((3, a) = B q,(A, (3) = (3
q,(a, (3)
=
0 elsewherewhere a and (3 refer to coordinates of poi~ts on t~e
boundary.
The Green's function for a rectangle is10
G(x, y, a, (3)
for
for
2 ~ . h
(m1rY). m1r ( )
- £-sm - - smh- b - (3
a m=l a a
. (m1rx) . (m1ra)
. sm ---;;- sm ~
m1r .
hm1rb
s m
-a a
2 ~.
(fn1r(3). m1r (b )
= - £-smh _.- smh- - y
a m=l a a
(23)
. (m1rx) . (m1ra)
. sm ~ sm
----;;-m1r .
hm1rb
s m
-a a
where
O:::;x:::;a
o :::;
{3 ~ b.The coordinates of the point being solved for are (x, y).
Equation 23 can also be expressed as a Fourier series in y by interchanging a and b, x and y, and a and {3. Both representations give the same results but, to obtain faster convergence the Fourier series in x should be used at points for which
(y: ~) > e ~ a)
(24)and the Fourier series in y when the reverse of Equation 24 is true. The exponential form of Equation 23 was used for developing the digital representations of the equation.
Using the Green's function for a rectangle and the two approximating regions 1234 and 4567 shown in Figure 8, the problem was programmed on an EAI 690 hybrid computer. The digital-analog problem split was as follows:
The digital computer calculated which of the two approximating regions would require the least random walks, the contribution the coincident boundary por-tions made to the solution, and the number of walks to be made from each non-coincident segment. In ad-dition, the digital computer controlled the modes of operation of the analog computer, monitored sense lines signaling boundary intersections, stored the boundary intersection coordinate values, computed the proper boundary function values, and set the initial conditions for the x and y integrators.
The analog computer performed the random walks by integrating the noise inputs, provided shaping filters for the noise, and with analog comparators and parallel
B 4 J
>t
\12
+ •
01
-I
Dr--~-
---
:r-1~---
20 (A.!D)/2 A
x
Figure 8-Region for which solution to Laplace's equation is desired
Variance Reduction Technique 25
TABLE I-Solution for Laplace's Equation for Region of Figure 8 with Base Walk Number of 1000
1000 Walks SOLUTION USING Ra WITHOUT Ra Y DIGITAL HYBRID IERRORI #WALKS HYBRID IERRORI
X = 100
36 35.33 35.19 .14 49 35.60 .28
32 30.62 30.54 .08 97 31.56 .94
28 25.83 2.5.89 .06 143 26.52 .31
26 23.38 23.55 .17 161 24.08 .70
24 20.89 20.72 .17 178 20.56 .33
22 18.33 18.25 .08 190 18.76 .43
20 15.69 15.63 .06 212 15.76 .07
18 12.93 13.20 .27 201 13.00 .07
16 10.02 9.89 .13 191 10.76 .75
14 6.90 7.11 .21 154 8.20 1.30
12 3.54 3.42 .12 91 4.48 .94
logic components simulated the boundary and flagged the digital computer when an intersection occurred.
Table I contains the results for the region shown in Figure 8 with dimensions:
A =200 B = 40 D = 10 H = 10,
a 5-unit boundary segment for region 1234, a 2-unit segment for region 4567, and a base walk number of
TABLE II -Solution for Laplace's Equation for Region of Figure 8 with Base Walk Number of 10,000
SOLUTION USING Ra
Y DIGITAL HYBRID I ERROR I # WALKS
x
= 10036 3.5.33 35.34 .01 490
32 30.62 30.61 .00 970
28 25.83 25.71 .12 143
26 23.38 23.40 .02 1610
24 20.89 20.82 .07 1780
22 18.33 18.33 .00 1900
20 15.69 15.66 .03 2120
18 12.93 12.99 .06 2010
16 10.02 10.02 .00 1910
14 6.90 6.93 .03 1540
12 3.54 3.61 .07 910
2€) Spring Joint Computer Conference, 1970
1000. The solutions both with and without an approxi-mating region are compared to the results obtained using finite differences and a relaxation type numerical solution.
Table II gives the results using an approximating region and a base walk number of 10,000. Solutions for many other points were madel l with similar reductions in error. Note in Tables I and II that less than 25% of the base number of walks is actually made.
QUICK LOOK CAPABILITY
One advantage of the continuous random walk with or without Ra is the ability to find the solution for just one point in the region while the solution by con-ventional digital techniques requires the solution over the whole region to obtain the solution at a particular point. Using the continuous random walk technique a solution for cJ> can be found at a point or points as the geometry is changed allowing a quick look at the effect of a geometry change at some critical point. For the problem of the preceeding section Figure 9 shows the resultant change in cJ> at the point x = 100, Y = 20, as H varies from 0 to 20 with a segment increment of 5 and 1000 for the· base number of walks. The data for the curve was taken in less than 15 minutes, the time required for one conventional digital solution.
20
1.5
+
10
.5
o .5 10 1.5 20
H
Figure 9-Solution of Laplace's equation for X = 100, Y = 20 in the region of Figure 8 as H is varied
SOLUTION OF POISSON'S EQUATION WITH CONSTANT FORCING FUNCTION
SOLUTION OF POISSON'S EQUATION WITH CONSTANT FORCING FUNCTION