Numerical integration: Trapezoidal method

We will calculate the area under the curve $f(x)$ within the bounds $x=a$ and $x=b$. We will divide the curve into $n$ trapezoids with width $h=(b-a)/n$. Area under the curve:

$$A_n = h \left[\frac{f(a)}{2} + f(a+h) + f(a+2h) + \cdots + \frac{f(b)}{2}\right]$$

src/10_trapezoid.f90

! Caculate the area of sine curve using trapezoidal method

PROGRAM trapezoid
  IMPLICIT none

  INTRINSIC :: sin
  REAL :: a, b, h, area
  INTEGER :: i, n

  n = 100
  a = 0
  b = 3.14

  h = (b - a) / n

  area = 0.5 * (sin(a) + sin(b))

  DO i = 1, n-1
    area = area + sin(a + i*h)
  END DO

  area = h * area

  PRINT *, "Area = ", area

END PROGRAM trapezoid

Numerical integration II

We will try to improve above numerical integration code. Previously we had the integration function hard coded, now we will pass it as argument from the calling program and use separate module for the calculation.

First we write our executable code in a module:

src/14_integration_module.f90

MODULE integration_module
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: integral

  CONTAINS

  FUNCTION integral(f, a, b, n) RESULT(integral_result)

    INTERFACE
      FUNCTION f(x) RESULT(f_result)
        REAL, INTENT(IN) :: x
        REAL :: f_result
      END FUNCTION f
    END INTERFACE

    REAL, INTENT(IN) :: a, b
    INTEGER, INTENT(IN) :: n
    REAL :: integral_result
    REAL :: h, area
    INTEGER :: i

    h = (b - a) / n

    area = 0.5 * (f(a) + f(b))

    DO i = 1, n-1
      area = area + f(a + i*h)
    END DO

    integral_result = h * area
  END FUNCTION integral
END MODULE integration_module

src/14_integration.f90

PROGRAM integration
  USE integration_module
  IMPLICIT NONE
  INTRINSIC :: sin

  PRINT *, integral(sin, a=0.0, b=3.14159, n=1000)
END PROGRAM integration