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: '
Print *, x0, x1
PRINT *

PRINT *, 'Number of panels (even): '
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
```

Come back to the previous page

Page builder: Charles Boivin