SIAM 2002 challenge # 10

This solution was supplied by Marijke van Gans in HTML form

First, let's anticipate some later results to get an idea of the magnitudes involved.

Above: vertical covers all the probability. Red too small to see.

Below: same plot, the 100 pixels vertical stretched to 1,000,000,000 pixels, and only showing the very bottom portion. Horizontal scale is the same as above.

Ten places to the right of the decimal point gives us six zeros and only four significant figures. Let's interpret Trefethen's "ten digits" the hard way and go for 10 significant digits, i.e. sixteen places. Preliminary integration at coarse Delta t steps shows that for that kind of precision we'll need to integrate up to t = 8 or 9, make it 10 for good measure.

Modelling finite random walks in 1 dimension

      (n C k) / 2ⁿ   where   (n C k) := n! / (k! (n-k)!)

Taking a step left or right makes 2s difference, so we are an integer number times 2s from (for even n) the origin, and that integer is n/2 - k. This binomial distribution has a mean of n/2 for k, i.e. a mean of 0 in terms of position. It has a standard deviation of sqrt(n)/2 for k, i.e. of sqrt(n)s in terms of position.

From a dimly remembered theorem: in the limit for large n, the binomial approaches a Gaussian (aka "normal") distribution. This too must then be the one with position mean at the origin and standard deviation sqrt(n)s.

Modelling infinitesimal random walks in 1 dimension

      std. dev.  =  sqrt(n) s  =  sqrt(wt) s

is kept finite, for all finite t. That means that, in the limit for s -> 0 and w -> oo, we must have w growing as 1/s². Let us assume that the limit of this behaviour is what is meant with "random walk with infinitesimal step lengths". It no longer resembles a Gaussian but is one. One consequence of ws now being infinite is that there is a small but nonzero probability of ending up anywhere. There's no longer a front of "furthest we could have gotten" at any one time.

For simplicity let us assume sqrt(w)s is not just finite in the limit but actually 1. This merely sets a scale of units on our t and makes no difference to the answer. Standard deviation is now sqrt(t).

Modelling infinitesimal random walks in 2 dimensions

Conversely, if a 2-D distribution is a Gaussian of the distance r to some point, then the density function can be written as the product of a function that only depends on x and one that only depends on y (Gaussians are the only function with this property).

Random walks with capture at the edges

Note that the creaming off also skews the remaining distribution inside! If any stuff that reaches "there" isn't allowed to play anymore then we see that reflected "here" as a biased thinning out. What exactly is the shape of the depletion? Well, in somewhat Huygensian fashion, each point of the "wave front" is again an unbiased source. Brownianly moving particles that reach the edge start, from there, a memory-less random walk.

This starts looking like some kind of functional equation.

We have the original "source" of stuff which, if unhindered, would diffuse to attain, at any time t, a density

      C_t exp(-r² / (2t))   =   C_t exp((-x²-y²) / (2t))
      where C_t   =   (2pi t)^-½

at the point x,y at time t. The real depleted density (summing to less than 1 over space) is some function that is depleted by an amount that depends on how much the function was everywhere earlier... Several strategies spring to mind. One is a hideously complicated spatio-temporal integral differential equation. Sadly, i'm not clever enough to do that kind of stuff. So let's try a different tack.

A Brownianly moving particle that reaches the edge has subsequently, if we didn't cream it off, equal probability drifting either side of the edge (further out or back in). So creaming off depletes both outside and inside the edge. In other words: some of the distribution we see inside, in an un-creamed off distribution, consists of Brownian movers that have been past the edge and drifted back. It is these returnees we must eliminate.

What is the shape of this depletion? Well, the tail outside the edge consists completely of movers that have crossed the edge. So that is cut off whole and we now see we must mirror that tail and subtract the same amounts also inside the edge (because of the 50/50 chance of drifting back).

So far so good... except that the two mirror tails eventually start to overlap noticeably. A more complete analysis along these lines leads to the following image (and we might as well state it in the 2-D case of the SIAM problem straightaway):

The solution

Is it the Right Thing to have non-physical sources two rectangles away that contribute in the positive??? Well yes, they constitute the depletion of the depletion , it all works out i think. [Someone better check whether my ramblings make any sense -- MvG].

Our probability density function has become

      W(t;x,y)   =   Sum_i Sum_j (-1)^{i + j} C_t² exp((-[x-i]²-[y-10j]²) / (2t))
      where C_t   =   (2pi t)^-½

which has all the depletion built in. Presumably it nicely evens out to zero over time. We still need to look at the depletion itself of course (in order to compare depletion at the short sides with depletion at the long sides) but at least we now know what the function at any time is, from which that depletion gets recruited.

Note that the density becomes zero when we reach the edges (because of symmetry with the negative rectangles). This looks counterintuitive and at first caused me to abandon the idea. But the quickie coarse grained simulation (yes, UBASIC again :) reproduced in the animated GIF above shows that this (density function graph intersecting zero at the edge) really happens!

This suggests we can't integrate the density at the edge over time to count the escapees along short and long edge (that density is zero). But here's a method that does work:

Note that our infinite-sum W can still be written as a product

      W(t;x,y)   =   X(t;x) Y(t;y)

where

      X(t;x)   =   Sum_i (-1)ⁱ C_t exp(-[x-i]² / (2t))

      Y(t;y)   =   Sum_j (-1)^j C_t exp(-[y-10j]² / (2t))

      where C_t   =   (2pi t)^-½

Now at each moment in time, the amount of probability still in the rectangle is HV where

      H(t)   =   Int[-½...+½] X(t;x) dx

      V(t)   =   Int[-5...+5] Y(t;y) dy

and the rate of escape over the long edge is

      dL/dt   :=   (-dH/dt) V(t)

and the rate of escape over the short edge is

      dS/dt   :=   H(t) (-dV/dt)

so we need the time integrals

      L(oo)   =   Int[0...oo] (-dH/dt) V dt

      S(oo)   =   Int[0...oo] H (-dV/dt) dt

Any lingering doubts whether this approach is correct evaporated when a quick cross check using a simulation in C yielded (to a few significant figures) the same result. This clinches it because that C program uses an entirely different approach (simply an N × 10N grid and finite steps in discrete time). Its results:

SPATIAL steps   TIME  PROBABILITIES
in 1  in 1/4    steps "blue" = prob. of  "red" = prob. of   "white" = prob.
unit  rectangle done  escape long side   escape short side  remaining
----  ------    ----  -----------------  -----------------  -----------------
   8   4x 40     600  0.999999538109659  0.000000461797641  0.000000000092700
  16   8x 80    2400  0.999999597667567  0.000000402227201  0.000000000105233
  32  16x160    9600  0.999999611579003  0.000000388312412  0.000000000108585
  64  32x320   38400  0.999999614997764  0.000000384892799  0.000000000109437
 128  64x640  153600  0.999999615848790  0.000000384041557  0.000000000109652
  oo extrapolation:                      0.0000003838

The program

For f = 1/sqrt(t), with change of variable u := xf, we have

      Int[a...b] (2pi t)^-½ exp(-[x-i]² / (2t)) dx   =   Int[af...bf] (2pi)^-½ exp(-[u-if]² / 2) du

in other words, we don't have to keep doing a new integral dx for each new t but it's the same old normal distribution each time, just the bounds keep shifting. For any one value of t, and hence of f, H is given by integrating (2pi)^-½exp(x²/2) from -½f to +½f and counting it positive, from -1½f to -½f and from +½f to +1½f and counting it negative, from -2½f to -1½f and from +1½f to +2½f and counting it positive again, and so on. For V, the cutoffs are ±5f, ±15f, and so on.

This means we can calculate dH, and dV, as just the difference between such sums i.e. it involves a really short piece of integration e.g. from ½f to ½(f + df). The total workload, for all t, is just one whole normal distribution integral (from 0 to infinity) for each index i or j. And given that (2pi)^-½exp(x²/2) peters out to nothingness for ±10 at our resolution, and our integral is all sewn up by t = 9, i.e. f = 1/3, we only need to look at 30 i-terms and 3 j-terms. Here's the program that implements that.

PS - yes i've noticed i'm initialising H, V, L and S twice <g>, i've left it in as this is what i actually submitted and because the first initialisation (for t = 0) carries the comment with explanation. The second (for t = 1/400, f = 20) is what we really need because starting  f at infinity would make the running time rather long...

    5   point 6
   16   Er=1/10^24
   63   O=pi(2)
   64   C=1/sqrt(O)
   91   H=1 ' Int[-0.5..+0.5] Horz(t;x) dx
   92   V=1 ' Int[-5...+5]    Vert(t;y) dy
   93       ' where prob density(t;x,y) = Horz(t;x) * Vert(t;y)
   94       ' so Int[]Int[] prob dens dx dy = H*V
   95   L=0 ' Int[0...T] (-dH/dt)*V dt (overflow on long side)
   96   S=0 ' Int[0...T] H*(-dV/dt) dt (overflow on short side)
  105   H = 0.999999999999999999999985
  106   V = 1
  107   L = 0.000000000000000000000015
  108   S = 0
  199      F=20  :  DF=1 '                     where F=1/sqrt(T)
  200   if F=2.0 then gosub 300 : DF=1/1024
  201   if F=1.5 then gosub 300 : DF=1/(2^20)
  202   if F=0.5 then gosub 300 : ' keep DF
  203   if F=0.3125 then gosub 300 : end
  210     Fold=F : F-=DF
  221     N=0:Nh=0.5:DHq=0 '    DHq is quarter DH
  222     repeat
  223       DG=fnSimpson(&fnGauss,Nh*F,Nh*Fold,Er)
  224       if odd(N) then DHq+=DG else DHq-=DG
  225       inc N : Nh += 1
  227     until DG < Er
  228     DHq+=DHq:DH=DHq+DHq ' DHq is half DH
  229     '
  231     N=0:N5=5:DVq=0 '      DVq is quarter dV
  232     repeat
  233       DG=fnSimpson(&fnGauss,N5*F,N5*Fold,Er)
  234       if odd(N) then DVq+=DG else DVq-=DG
  235       inc N : N5 += 10
  237     until DG < Er
  238     DVq+=DVq:DV=DVq+DVq ' DVq is half DV
  239     '
  240     L-=DH*(V+DVq) ' Delta H * avg V
  241     S-=(H+DHq)*DV ' avg H * Delta V
  242     H+=DH
  243     V+=DV
  290     goto 200
  299   '
  300   print using(2,3),1/(F*F);
  301   print using(2,18),L;S;H*V;
  302   print using(3,1),Nh;
  303   print using(3,1),N5
  304   return
  305   '
  750   fnGauss(Z)
  760   return(C*exp(-Z*Z/2))
  799   '
  800   'Simpson's Rule (Numerical Integration)
  801   'adapted from SIMPSON.UB bundled with Yuji Kida's UBASIC,
  802   'after S.Moriguti Suuti Keisanjutu (Kyouritu Shuppan) p199
  803   '
  810   fnSimpson(&fnF(),A,B,Err)
  811      local D1,D2,Arg,M,N,S1,S2,I,I_I
  820   N=1:D1=B-A
  830   S2=(fnF(A)+fnF(B))/2:S1=0
  840   I=S2*2*D1/3
  850   repeat
  860      N*=2:D2=D1:D1/=2:S2+=S1:S1=0
  870      Arg=A+D1:for M=1 to N step 2:S1+=fnF(Arg):Arg+=D2:next
  880      I_I=I:I=(2*S1+S2)*D2/3:I_I-=I
  890   until absmax(I_I)<Err
  899   return(I)

--marijke

<back