Anonymous user
Numerical integration: Difference between revisions
→{{header|Fortran}}
Line 498:
=={{header|Fortran}}==
In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:
ELEMENTAL FUNCTION ELEMF(X)
REAL :: ELEMF, X
ELEMF = F(X)
END FUNCTION ELEMF
Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function
PROGRAM INTEGRATE
INTEGER, PARAMETER :: N = 20 ! or whatever
REAL, PARAMETER :: A = 0.0, B = 15.0 ! or whatever
REAL, PARAMETER :: H = (B - A) / N
REAL, PARAMETER, DIMENSION(0:2*N) :: XPOINTS = (/ (A + H*i/(N*2), i = 0, 2*N) /)
REAL, DIMENSION(0:2*N), TARGET :: FPOINTS ! gather up all the f(x) values needed for all methods once and distribute via pointers
REAL, DIMENSION(:), POINTER :: FLEFT, FMID, FRIGHT
REAL :: LEFTRECT, MIDRECT, RIGHTRECT, TRAPEZOID, SIMPSON
FPOINTS = ELEMF(XPOINTS) ! elemental function wrapper for F runs once for each element of XPOINTS
FLEFT => FPOINTS(0 : 2*N-2 : 2)
FMID => FPOINTS(1 : 2*N-1 : 2)
FRIGHT => FPOINTS(2 : 2*N : 2)
! Left rectangular rule integral
LEFTRECT = H * SUM(FLEFT)
! Middle rectangular rule integral
MIDRECT = H * SUM(FMID)
! Right rectangular rule integral
RIGHTRECT = H * SUM(FRIGHT)
! Trapezoid rule integral
TRAPEZOID = H / 2 * SUM(FLEFT + FRIGHT)
! Simpson's rule integral (coming soon)...
SIMPSON = H / 6 * SUM(FLEFT + FRIGHT + 4*FMID)
END PROGRAM INTEGRATE
=={{header|Haskell}}==
| |||