program test2dintegrate implicit none double precision :: x,y,r,theta,pi double precision :: xmax,ymax,rmax,rl,rh,A,d,IMC,Ireference,NMC,Nrepeat integer :: i,k,l,ix,iy,n,ninterval,nininterval,npoints double precision, external :: stick double precision, external :: grnd pi=2.0d0*asin(1.0d0); ! Integration interval in cartesian coordinates xmax=3.0d0; ymax=3.0d0; ! Integration interval in polar coordinates rmax=3.0d0; ! Initialize MC runs call sgrnd(1782) NMC=900; Nrepeat=10; ! Part b): Direct MC in cartesian do n=1,Nrepeat IMC=0.0d0 do i=1,NMC x=(2.0d0*grnd()-1.0)*xmax; y=(2.0d0*grnd()-1.0)*ymax; r=sqrt(x*x+y*y); theta=atan2(y,x); IMC=IMC+stick(r,theta) enddo print *,'Direct xy MC estimate of integral is',2*xmax*2*ymax*IMC/NMC enddo ! Part c): Direct MC in polar coordinates do n=1,Nrepeat IMC=0.0d0 do i=1,NMC ! Should weight correctly with area element r dr dtheta ! i.e. generate points distributed as r ! Integral of r from 0 to rmax is 1/2 rmax^2 ! This gives normalized function r/(1/2 rmax^2). ! and hence normalized integrated function r^2/rmax^2 ! Solving u=r^2/rmax^2 gives r=sqrt(u*rmax^2) r=sqrt(rmax**2*grnd()); theta=2*pi*grnd(); IMC=IMC+stick(r,theta) enddo print *,'Direct polar MC estimate of integral is',pi*rmax**2*IMC/NMC enddo ! part d): fully stratified MC in cartesian coordinates do n=1,Nrepeat npoints=0 ninterval=int(sqrt(NMC)) IMC=0.0d0 do k=0,ninterval-1 do l=0,ninterval-1 x=(k-ninterval/2+grnd())/(ninterval/2)*xmax; y=(l-ninterval/2+grnd())/(ninterval/2)*ymax; r=sqrt(x*x+y*y); theta=atan2(y,x); IMC=IMC+stick(r,theta) npoints=npoints+1 enddo enddo print *,'Fully stratified xy MC estimate of integral is',2*xmax*2*ymax*IMC/npoints enddo ! part e): partially stratified MC in polar coordinates do n=1,Nrepeat ninterval=20; npoints=0; IMC=0.0d0 do k=0,ninterval-1 ! r interval to handle rl=k*rmax/ninterval; rh=(k+1)*rmax/ninterval; ! Figure out number of points needed in this interval ! note that because of roundoff errors final total number of points ! may not be exactly = NMC A=pi*rh**2-pi*rl**2 nininterval=nint(A/(pi*rmax**2)*NMC) do l=1,nininterval ! Should weight correctly with area element r dr dtheta ! i.e. generate points distributed as r ! Integral of r from rl to rh is 1/2 (rh^2 - rl^2). ! This gives normalized function r/(1/2(rh^2 - rl^2)). ! and hence normalized integrated function (r^2-rl^2)/(rh^2 - rl^2) ! Solving u=(r^2-rl^2)/(rh^2 - rl^2) gives r=sqrt(u*(rh^2 - rl^2)+rl^2) r=sqrt(grnd()*(rh**2 - rl**2)+rl**2); ! Following would be WRONG: !r=grnd()*(rh-rl)+rl; theta=2*pi*grnd(); IMC=IMC+stick(r,theta) npoints=npoints+1 enddo enddo print *,'Stratified polar MC estimate of integral is',pi*rmax**2*IMC/npoints enddo ! Part a): get accurate answer with very high precision ! numerical integration Ireference=0.0d0; ! Loop over increasingly accurate d values d=0.1d0 do npoints=0 ! Note the deliberate 01 at the end for comparison roundoff avoidance! if (d < 0.000101) exit n=nint(xmax/d) do ix=-n,n do iy=-n,n x=ix*d; y=iy*d; r=sqrt(x*x+y*y); theta=atan2(y,x); Ireference=Ireference+stick(r,theta) npoints=npoints+1 enddo enddo ! Ireference=Ireference*d*d ! Correspondence to MC integration is clearer if using: ! Note the use of n+0.5; draw a grid and you'll understand why Ireference=((n+0.5d0)*d*2)**2*Ireference/npoints print *,'Using d=',d,' reference value is',Ireference d=d/10.0d0; enddo end program test2dintegrate double precision function stick(r,theta) implicit none double precision :: r,theta,pi,f,d pi=2.0d0*asin(1.0d0); d=3.0d0/sqrt(2.0d0)+0.5d0*cos(theta*3)**2; if (r < d-0.5d0) then f=1.0; else if (r < d+0.5d0) then f=0.5d0-0.5d0*sin(pi*(r-d)) else f=0.0d0; endif stick=f end function stick !************************************************************************ ! Mersenne twister random number generator !************************************************************************ subroutine sgrnd(seed) ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) ! dimension mt(0:N-1) ! the array for the state particle common /block/mti,mt save /block/ ! ! setting initial seeds to mt[N] using ! the generator Line 25 of Table 1 in ! [KNUTH 1981, The Art of Computer Programming ! Vol. 2 (2nd Ed.), pp102] ! mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue ! return end !************************************************************************ double precision function grnd() ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) ! constant vector a parameter(UMASK = -2147483648) ! most significant w-r bits parameter(LMASK = 2147483647) ! least significant r bits ! Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) ! dimension mt(0:N-1) ! the array for the state vector common /block/mti,mt save /block/ data mti/N1/ ! mti==N+1 means mt[N] is not initialized ! dimension mag01(0:1) data mag01/0, MATA/ save mag01 ! mag01(x) = x * MATA for x=0,1 ! TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) ! if(mti.ge.N) then ! generate N words at one time if(mti.eq.N+1) then ! if sgrnd() has not been called, call sgrnd(4357) ! a default initial seed is used endif ! do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif ! y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) ! if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif ! return end