### MC-INTEGRATION2

! Program name: MC-INTEGRATION2

! INTEGRATION BY MONTE CARLO METHOD

! INTEGRAL BY MEAN FUNCTION

! EXAMPLE FUNCTION f(x,y)=4-x^2-y^2

! In the interval (0,5/4) and (0,5/4) the integral is 4.62239

Program monte_carlo

real :: f, f2, integrand, integral, sigma, a, b, x, c, d, y

integer :: i, n

data exact/4.62239/

Print *, "Enter the grid of the integration, n: "

If (n .eq. 0) stop

Print *, "Enter the lower limit and upper limit of first integration, a, b:"

Read *, a, b  ! In the suggested problem a=2 and b=5/4

Print *, "Enter the lower limit and upper limit of second integration, c, d:"

Read *, c, d  ! In the suggested problem a=0 and d=5/4

f= 0.0

f2=0.0

call random_seed()

do i=1,n

call random_number(x)

x=x*b

call random_number(y)

y=y*d

f= f+integrand(x,y)

f2=f2+(integrand(x,y)**2)

!   print *, x, f, f2

End do

f=f/n

f2=f2/n

integral=(b-a)*(d-c)*f

sigma=(b-a)*(d-c)*sqrt((f2-f**2)/n)

Print *,

Print *, integral, sigma, "error= ", exact-integral

End program monte_carlo

function integrand(x,y) result(func)

Implicit none

Real :: func, x, y

func= 4-x**2-y**2

End function integrand