Numerical Integration

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)

\begin{align} I=\int_0^3 {{x^4 e^x} \over{ (e^x - 1)^2 }} {\mathrm d}x \end{align}

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

is mostly stuff we've seen before. The goal is to print out a table of x and f(x) for us to look at. However, the write statement inside the do loop is now calling a 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

tells the compiler that we will begin defining 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

all of this should look familiar to you—-it is just a block of F90 code that is labeled as a function. It takes 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

tells the compiler that this is the end of the program. Note that the name here "trapezoid0" has to match the name defined at the beginning of the program. Don't worry, the compiler will tell you!

---

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

What this does is say 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:

1. 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.
2. 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)

First, we 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)

A 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