P111.F90

Demo of various integration techniques

!
! =====> Program - P124.F90
!
        PROGRAM Intgrt
!  Program to demonstrate numerical integration.


!       Define function return types!!
        REAL            MidPt, Trapzd, Simpsn

        INTEGER         Panels
        REAL            x0, x1

        PRINT *, 'This is Program >> P124 = Integration Demo'
!
!     Tell program where data for  READ *  is coming from
      OPEN(UNIT=5, FILE='P124.DAT')      ! UNIT=5 is the default input
!
        PRINT *, 'Demonstration of Numerical Integration'
        PRINT *,'Function = 1/(1 + x*x)  for example'
        PRINT *
        PRINT *, 'Low limit, high limit: '
        READ *, x0, x1
        Print *, x0, x1
        PRINT *

        PRINT *, 'Number of panels (even): '
        READ *, Panels
        Print *, Panels
        PRINT *

        PRINT *, 'Midpoint:  ', MidPt(x0, x1, Panels)
        PRINT *, 'Trapezoid: ', Trapzd(x0, x1, Panels)
        PRINT *, 'Simpson''s  ', Simpsn(x0, x1, Panels)

        STOP
        END


        REAL FUNCTION F(x)
!       The function to be integrated.

        F = 1/(1 + x*x)
        RETURN
        END


        REAL FUNCTION MidPt(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using the Midpoint Rule with N panels.

        REAL            x0, x1,    &
                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = 0.0

!       Add area of each panel, assuming unit width
!       and multiplying by the true width h later.
        x = x0+h/2.0   ! Initial start is middle of first panel
L1:     DO k= 1,N
          Sum = Sum + F(x)
          x = x + h
        END DO L1

        MidPt = h*Sum

        RETURN
        END


        REAL FUNCTION Trapzd(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using the Trapezoidal Rule with N panels.

        REAL            x0, x1,    &
                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = (F(x0) + F(x1))/2.0

!       Add contributions from interior trapezoids.
        x = x0+h
L1:     DO k = 1,N-1
          Sum = Sum + F(x)
          x = x + h
        END DO L1

        Trapzd = h*Sum

        RETURN
        END


        REAL FUNCTION Simpsn(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using Simpson's Rule with N/2 panels.

        REAL            x0, x1,    &
                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = F(x0) + F(x1)

!       Add 4 * odd-indexed f values and 2 * even-
!       indexed f values to Sum
        x = x0
L1:     DO I = 1, N-1
          x = x + h

          IF (Mod(I, 2) == 0) THEN
            Sum = Sum + 2.0*F(x)
          ELSE
            Sum = Sum + 4.0*F(x)
          END IF
        END DO L1

        Simpsn = h/3.0*Sum

        RETURN
        END

DATA:
-2.0  2.0
10

OUTPUT:
Program entered
 This is Program >> P124 = Integration Demo
 Demonstration of Numerical Integration
 Function = 1/(1 + x*x)  for example

 Low limit, high limit: 
  -2.0000000   2.0000000

 Number of panels (even): 
 10

 Midpoint:     2.2164154
 Trapezoid:    2.2100482
 Simpson's     2.2150481
Fortran-90 STOP

Page builder: Charles Boivin

Last modified: 11/07/95