! ! =====> 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