The following solution was supplied by Brian Medley. In private conversation he remarks that he only coded it up as a finite simulation (rather than taking the analytical route) in order to provide a cross check, as there was an answer on the board already. Nevertheless, this approach also yields more than ten correct digits.
SIAM Problem 6.
Essentially, I set p(start)=1 and p(everywhere else)=0 in a big array, then update all the probabilities for each timestep. After every step, I reset p(start) to 0 (so that a flea that returns plays no further part in the process). I keep track of: p=cumulative probability of returning to the start so far, pleft=total of all the "live" probabilities still sloshing around the array and 3 other quantities that tot up how much probability "leaks" over each edge of the array. The total of all 5 quantities is 1.
I found that with a large enough array, practically no probability escapes to the West ("upwind"), which justifies the assumption that the probability that escapes to the East is never coming back (East and West boundaries are equidistant from the start). And we can get a small enough probability disappearing North and South rather than East to reason that by the time any of this probability comes back again - without disappearing over the East border first - it, too, will have dwindled to practically zero.
_________________________________________________________________________ 10 'global initialisation 20 ' 30 point 7 40 emaword 7 50 ' 60 ' Imax and Jmax set the dimensions of an array 70 ' of points for which we will track the probability 80 ' of the photon being present at t=0, 1, 2,... 81 ' Imax=64 and Jmax=100 seem high enough to get Ep to 10 s.f. 90 ' 100 Imax=26:'must be an even number 110 Jmax=40 120 dim ema(1;2*Imax,Jmax) 130 ' 140 ' The next two lines set Ep0 and Ep1 which are two values we 150 ' happen to know epsilon is between. Initially they could be 0 and 1/4 160 ' 170 Ep1=0.25 180 Ep0=0.0 190 ' 200 '---------------START WITH A NEW Ep----------------- 210 ' 220 ' We hone in on epsilon by interval bisection, taking a new Ep 230 ' halfway between Ep0 and Ep1. 240 ' 250 Ep=(Ep1+Ep0)/2 260 ' 270 ' The variables on the next line are used to track how probability 280 ' "leaks" out of the array, helping estimate error. 290 ' 300 Eleak=0:Wleak=0:Nsleak=0 305 print 310 print "trying";Ep 320 ' 330 ' The next six lines initialise the array. The flea starts at Imax,0 340 ' not 0,0 to avoid negative array subscripts. We use symmetry to avoid 350 ' negative j's: P(i,j)=P(i,-j) so no need to store both. 360 ' 370 clr block ema(1;*,*) 380 ' 390 for I=0 to 2*Imax 400 for J=0 to Jmax 410 ema(1;I,J)=0 420 next J 430 'next I 440 ' 450 ema(1;Imax,0)=1:'starting position of the flea. 460 ' 470 P=0:'used to store p(flea has returned to start) 480 N=0:'number of stages or time steps. I could have called this t. 490 ' 500 '--------------START OF MAIN LOOP-------------- 510 ' 520 Pleft=0:'used to see how much probability is still sloshing around 530 F=N@2:'F records whether N is even or odd 540 ' 550 ' F is important as I use the fact that at any time step, half the 560 ' probabilities in the array (in a chessboard pattern) are zero. 570 ' F tells us which half: at each stage we overwrite the old zeros 580 ' with new proabilities: later we will overwrite the old non-zero 590 ' probabilities with zero. 600 ' 610 ' The next nine lines update probabilities in the row j=0 620 ' This row is special as ema(1;K,1) is doing double duty for 630 ' positions (K,1) and (K,-1), both of which contribute to the 640 ' next probability for (K,0) 650 ' 660 for K=0 to 2*Imax 670 if K@2=F then goto 740 680 Pr=0 690 if K>0 then Pr=Pr+ema(1;K-1,0)*(1/4+Ep) 700 if K<2*Imax then Pr=Pr+ema(1;K+1,0)*(1/4-Ep) 710 Pr=Pr+ema(1;K,1)/2:'/2 not /4 because (K,1) is doing double duty 720 Pleft=Pleft+Pr 730 ema(1;K,0)=Pr 740 next K 750 ' 760 ' now update probabilities in the rest of the array 770 ' 780 for J=1 to Jmax 790 for I=0 to 2*Imax 800 Pr=0 810 if (I+J)@2=F then goto 880 820 if I>0 then Pr=Pr+ema(1;I-1,J)*(1/4+Ep) 830 if I<2*Imax then Pr=Pr+ema(1;I+1,J)*(1/4-Ep) 840 Pr=Pr+ema(1;I,J-1)/4 850 if J<Jmax then Pr=Pr+ema(1;I,J+1)/4 860 Pleft=Pleft+Pr*2 870 ema(1;I,J)=Pr 880 next I 890 next J 900 ' 910 ' Finally, set all the old probabilities to zero. 920 ' 930 for J=0 to Jmax 940 for I=0 to 2*Imax 950 if (I+J)@2=F then ema(1;I,J)=0 960 next I 970 next J 980 ' 990 ' Update p and print it out 1000 ' 1010 P=P+ema(1;Imax,0) 1020 'print N;:print using(1,14),P,Pleft,Ep 1030 ' 1040 ' Check the result. If P rises above 0.5, Ep was too small 1050 ' whilst if P+Pleft is below 0.5, p can never reach 0.5 and so 1060 ' Ep was too big. 1070 ' 1080 if P>0.5 then Ep0=Ep:print "Ep was too small":gosub 1630:goto 200 1090 if P+Pleft<0.5 then Ep1=Ep:print "Ep was too big":gosub 1630:goto 200 1100 ' 1110 ' Reset the probability for the starting position to 0. 1120 ' This prob has been included in P and we are not interested in the 1130 ' flea's subsequent activities. 1140 ' 1150 ema(1;Imax,0)=0 1160 N=N+1 1170 'gosub 1000: 'a routine to show the current probs for a few cells 1180 gosub 1430:'a routine to find out how much probability is leaking away 1190 ' 1200 ' Fleas blown off the east boundary are assumed not to come back 1210 ' Other leaks are more serious and suggest errors if too high 1220 ' 1230 Eleak+=East:Wleak+=West:Nsleak+=Ns 1240 'print using(1,15),Eleak;Wleak;Nsleak 1250 goto 500 1260 ' 1270 '----------------END OF MAIN LOOP--------------- 1280 ' 1290 ' 1300 ' 1310 ' 1320 ' 1330 ' Subroutine to print out some probabilities near the starting point 1340 ' 1350 for J=0 to 4 1360 for I=-4 to 4 1370 print using(1,5),ema(1;I+Imax,J); 1380 next I 1390 print 1400 next J 1410 A$=input$(1) 1420 return 1430 ' 1440 ' Subroutine to calculate how much probability is leaking away over the 1450 ' edges of the grid. 1460 ' 1470 East=0:West=0:Ns=0 1480 ' 1490 for J=0 to Jmax 1500 East=East+ema(1;2*Imax,J) 1510 West=West+ema(1;0,J) 1520 next J 1530 ' 1540 East=(East*2-ema(1;2*Imax,0))*(1/4+Ep) 1550 West=(West*2-ema(1;0,0))*(1/4-Ep) 1560 ' 1570 for I=0 to 2*Imax 1580 Ns=Ns+ema(1;I,Jmax) 1590 next I 1600 ' 1610 Ns=Ns/2 1620 return 1630 ' 1640 ' print out values for last Ep 1650 ' 1660 print using(1,18),"Probablity of flea returning to start >=";P 1665 print using(1,18),"Probability of";Pleft;" still in array" 1670 print "Leakages:" 1671 print "East edge (normal exit, blown downwind)";Eleak 1672 print "West edge (exit upwind: this stuff will come back!)";:print using(1,18),Wleak 1673 print "N/S edges (drifting across the wind)";Nsleak 1680 return