In this chapter, we will try to solve **Example G-2** from MathChapter G ("Numerical Methods") from McQuarrie & Simon. In addition to showing an example of the trapezoidal and Simpson's rule solutions. To do this, we will learn two new features of F90: **suborutines/functions** and **if/then logical statements**

Recall (p. 633 of McQ&S) that we are trying to solve the integral:

(1)For simplicity, I will only do the n=10 solution (comparing to the table in Example G-2), and leave the additional calculations as an exercise to the reader

Before we begin, it is always a good idea to take a look at the function by computing the values that we will use. This is a way for you to make sure that you have implemented the function correctly. We can do this with a few lines of code

```
program trapezoid0
implicit none
integer, parameter :: n=10
real :: u
integer :: i
do i=0,n
u = 3.0*i/n
write (*,*) u, integrand(u)
end do
STOP
contains
function integrand(x) result (value)
implicit none
real :: x
real :: value
value = (x**4)*EXP(X)/((EXP(X)-1.0)**2)
end function integrand
end program trapezoid0
```

Let's take this in sections.

```
program trapezoid0
implicit none
integer, parameter :: n=10
real :: u
integer :: i
do i=0,n
u = 3.0*i/n
write (*,*) u, integrand(u)
end do
STOP
```

**function**named "integrand" to compute the value of the integrand at each point. We will print out a

This is specified in the section below:

`contains`

**functions**(pieces of code that return a number) and

**subroutines**(pieces of code that do some operation, but

*don't*return a number). These are only called if the code above "contains" calls them. Otherwise, they sit waiting patiently to be called upon. You can define as many functions and subroutines as you like. Here we have defined one function:

```
function integrand(x) result (value)
implicit none
real :: x
real :: value
value = (x**4)*EXP(X)/((EXP(X)-1.0)**2)
end function integrand
```

*x*as an input, and returns

*value*. You can of course name these whatever you like. Be sure to define them (as we have done here, we have made them both real) Note that the name of the argument for the function (here it is

*x*, above it was

*u*, do

**not**have to be the same.

**EXP**is simply a built-in F90 function to compute $e^x$.

Finally,

`end program trapezoid0`

Now, let's compile and run this program. Here's how I did it:

```
jschrier@ws7:~/pchem1/examples$ pgf90 -o trapezoid0.r trapezoid0.f90
jschrier@ws7:~/pchem1/examples$ ./trapezoid0.r
0.000000 NaN
0.3000000 8.9328051E-02
0.6000000 0.3493916
0.9000000 0.7574701
1.200000 1.278965
1.500000 1.871659
1.800000 2.490566
2.100000 3.092584
2.400000 3.640332
2.700000 4.104763
3.000000 4.466421
```

**Oops!**"NaN" is the compiler's way of telling us that we have divided-by-zero. The compiler doesn't understand limits. (You might demonstrate to yourself, by application of L'Hopital's rule that this should go to zero as x approaches zero.) We can get around this by modifying the function:

```
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
```

**if**"x" is less-than (

**.lt.**) $10^{-5}$,

**then**set x= $10^{-5}$. In other words, to avoid the division by zero, we can set x to some small epsilon. You might play around to see how the value of this epsilon (here $ 10^{-5} $]]) makes a difference. This brings up two points:

- Within an
**if (… ) then ….. end if**block, we can have as many lines as wel like (e.g., here we just have a single line). These lines of code can call other subroutines/functions, contain do-loops, etc. *Comparisons*: We used ".lt.". But of course other numerical comparisons can be used. These are:**.eq.**,**.ne.**,**.lt.**,**.le.**,**.gt.**,**.ge.**. What do you think they mean?

Using the code above, I changed the name of the program to trapezoid1, and compiled/ran it as follows:

```
jschrier@ws7:~/pchem1/examples$ pgf90 -o trapezoid1.r trapezoid1.f90
jschrier@ws7:~/pchem1/examples$ ./trapezoid1.r
0.000000 9.9729933E-11
0.3000000 8.9328051E-02
0.6000000 0.3493916
0.9000000 0.7574701
1.200000 1.278965
1.500000 1.871659
1.800000 2.490566
2.100000 3.092584
2.400000 3.640332
2.700000 4.104763
3.000000 4.466421
```

*Mission Accomplished!* We have fixed the divide-by-zero problem.

#### Trapezoid rule integration:

With a little modification, we can use this to evaluate the integral using the trapezoid rule, described on p.631 in McQuarrie and Simon. All we have done is insert a line to call a subroutine, and specified what that subroutine does.

```
program trapezoid2
implicit none
integer, parameter :: n=10
real :: u
integer :: i
do i=0,n
u = 3.0*i/n
write (*,*) u, integrand(u)
end do
call trapezoid_integration(n,3.0)
contains
subroutine trapezoid_integration(n,end_val)
implicit none
integer :: n
real :: end_val
real :: integral,u,h
integer :: i
integral = 0.0
do i=0,n
u = (end_val*i)/n
! implement Eqn G.4
if ((i.eq.0).or.(i.eq.n)) then
integral = integral+integrand(u)
else
integral = integral+(2.0*integrand(u))
end if
end do
h=end_val/n
integral = (h/2.0)*integral
write (*,*) '#trapezoidal integration = ',integral
end subroutine trapezoid_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 trapezoid2
```

What have we done?

`call trapezoid_integration(n,3.0)`

**call**a

**subroutine**named

*trapezoid_integration(n,end_val)*to actually perform the integration. The arguments are the number of points we are using to divide our function (

*n*) and the upper limit of integration (

*end_val*); the lower limit of integration is fixed at zero.

`subroutine trapezoid_integration(n,end_val)`

**subroutine**is like a function, but doesn't return a value. So there is not a need to specify a

**result**, as with the functions as described above. What is actually in this subroutine?

```
subroutine trapezoid_integration(n,end_val)
implicit none
integer :: n
real :: end_val
real :: integral,u,h
integer :: i
integral = 0.0
do i=0,n
u = (end_val*i)/n
! implement Eqn G.4
if ((i.eq.0).or.(i.eq.n)) then
integral = integral+integrand(u)
else
integral = integral+(2.0*integrand(u))
end if
end do
h=end_val/n
integral = (h/2.0)*integral
write (*,*) '#trapezoidal integration = ',integral
end subroutine trapezoid_integration
```

If the lower limit of the integral isn't 0:

-Need to define start_val (along with end_val). When the starting value is 0, we don't need to define it.

-Need to define u as (i*(end_val-start_val))/n + start_val

-Need to define h as (end_val-start_val)/n

#### Simpson's rule integration

From the discussion above, it should now be obvious to you how to implement Simpson's rule, as discussed on p.632 of McQuarrie & Simon.

Finally, I'll show you another approach which is not in the textbook…Monte Carlo integration