Here we will also introduce Monte Carlo integration.

Recall from basic calculus the mean value theorem for integrals: For $f(x)$ continuous over $[a,b]$, there exists a number, $c$, with $a < c < b$, such that

(1)Another way to think of this is that the area under the curve is the base width $(b-a)$ times the "average height" $f(c)$, in other words, we can rearrange the above formula into

(2)In some sense, the trapezoidal integration we studied in the previous section is nothing other than a set of midpoint integrations. (Think about it for a minute). Another approach we can take is to try to statistically estimate the midpoint by choosing points at random. In one dimension, this is actually kind of slow and innaccurate. But in higher dimensions it becomes one of the few practical methods of determining integrals.

To state the algorithm concisely, we will:

- Pick $n$ randomly distributed points $x_1, x_2, x_3, ...x_n$ in the interval of $[a,b]$.
- Determine the average value of the function $\langle f \rangle = {\frac{1}{n}} \sum_{i=1}^{n} f (x_i)$
- Compute the approximation to the integral using the mean value theorem $\int _a^b f(x)\,{\mathrm d}x \approx (b-a) \times \langle f \rangle$
- As a consequence of the Law of Large Numbers and/or the Central limit theorem, we can also estimate the error inherent in this approach, as $\mathrm{Error} \approx (b-a) {\sqrt{ {{ \langle f^2 \rangle - {\langle f \rangle }^2}\over {n}}}$, where $\langle f^2 \rangle = {\frac{1}{n}} \sum_{i=1}^{n} f^2 (x_i)$.

Let's consider again the case we studied in the previous section,

(3)we will use the same function for the integrand as before.

```
program mc_demo
implicit none
integer, parameter :: n=10000
call random_seed() ! needed to initialize the random number generator used in MC eval
call MC_integration(n,3.0)
contains
subroutine MC_integration(n,end_val)
implicit none
integer, intent :: n
real, intent :: end_val
real :: x, integral, integral_err
real (kind=8) :: f, f2
integer :: i
integral = 0.0
f= 0.0d0
f2 = 0.0d0
do i=1,n
call random_number(x)
x=x*end_val ! random_number only returns uniformly distributed from [0.0, 1.0]
f = f+integrand(x)
f2 = f2+ (integrand(x)**2)
end do
f=f/n
f2=f2/n
integral = (end_val-0.0)*f
integral_err = (end_val-0.0)*SQRT((f2 - f**2)/n)
write (*,*) "# MC integration = ",integral,"+/-",integral_err
end subroutine MC_integration
function integrand(x) result (value)
implicit none
real :: x
real :: value
if (x .lt. 0.00001) then
x = 0.00001
end if
value = (x**4)*EXP(X)/((EXP(X)-1.0)**2)
end function integrand
end program mc_demo
```

Let's take a look at the new parts:

```
integer, parameter :: n=10000
call random_seed() ! needed to initialize the random number generator used in MC eval
call MC_integration(n,3.0)
```

*random number generator*by calling one of F90's built-in subroutines "

**call random_seed()**". This built in random number (in reality, a pseudorandom number) generator is tolerable for the purposes of the class…better ones have been designed by computer scientists.

We then call the subroutine *MC_integration*. We first define a few variables, as you've seen before. One twist is that we want to keep extra bits of precision when we compute $f$ and $f^2$. To do this we can define these as

`real (kind=8) :: f, f2`

After defining the variables, the heart of this subroutine is:

```
integral = 0.0
f= 0.0d0
f2 = 0.0d0
do i=1,n
call random_number(x)
x=x*end_val ! random_number only returns uniformly distributed from [0.0, 1.0]
f = f+integrand(x)
f2 = f2+ (integrand(x)**2)
end do
```

- Get a random number using the F90 built in subroutine "
**call random_number(x)**". This puts a uniformly distributed number from $[0, 1]$ in the variable x. - We then have to multiply this to scale it over the base width of our integral
- We then keep running sums of $f$ and $f^2$.

Finally, after repeating this for $n$ times, we evaluate our averages, and write out the value for the integral and its error:

```
f=f/n
f2=f2/n
integral = (end_val-0.0)*f
integral_err = (end_val-0.0)*SQRT((f2 - f**2)/n)
write (*,*) "# MC integration = ",integral,"+/-",integral_err
```