Fortran numerical integration and root finding
last modified November 14, 2026
Two problems dominate practical numerical computing: computing an integral and finding where a function is zero. Fortran has no built-in routine for either, but both are short loops once the mathematics is clear, and both are good exercises in passing a function to another routine. Quadrature is the traditional name for numerical integration, and a root of a function is an argument at which the function returns zero.
This article implements the trapezoid, midpoint and Simpson rules and compares their accuracy on a table of interval counts, refines an integral until the answer stops changing, finds a root by bisection, by Newton's method and by the secant method, and finishes with Gaussian quadrature, which reaches the accuracy of a thousand intervals with five evaluations. Every rule receives the function it integrates as a procedure argument, with an explicit interface.
The trapezoid rule
The trapezoid rule replaces the area under the curve by a row of trapezia: the
interval from a to b is cut into n equal
parts of width h, and each part is approximated by a straight line
between the values of the function at its two ends. Summing the areas gives
h * (f(a)/2 + f(a+h) + ... + f(b)/2), where the two end values
count once and every inner value counts twice, because each is shared by two
neighbouring trapezia. The module below holds the rule together with its
relatives; its dummy argument is described by an abstract interface, so
the compiler checks that the function passed in has the right signature.
module rules
use iso_fortran_env, only: real64
implicit none
private
public :: trapezoid, midpoint, simpson, gauss2, gauss5
abstract interface
pure real(real64) function integrand(x) result(y)
import :: real64
real(real64), intent(in) :: x
end function integrand
end interface
contains
pure real(real64) function trapezoid(f, a, b, n) result(s)
procedure(integrand) :: f
real(real64), intent(in) :: a, b
integer, intent(in) :: n
real(real64) :: h
integer :: i
h = (b - a) / real(n, real64)
s = 0.5_real64 * (f(a) + f(b))
do i = 1, n - 1
s = s + f(a + real(i, real64) * h)
end do
s = s * h
end function trapezoid
pure real(real64) function midpoint(f, a, b, n) result(s)
procedure(integrand) :: f
real(real64), intent(in) :: a, b
integer, intent(in) :: n
real(real64) :: h
integer :: i
h = (b - a) / real(n, real64)
s = 0.0_real64
do i = 1, n
s = s + f(a + (real(i, real64) - 0.5_real64) * h)
end do
s = s * h
end function midpoint
pure real(real64) function simpson(f, a, b, n) result(s)
procedure(integrand) :: f
real(real64), intent(in) :: a, b
integer, intent(in) :: n
real(real64) :: h
integer :: i
h = (b - a) / real(n, real64)
s = f(a) + f(b)
do i = 1, n - 1
if (mod(i, 2) == 0) then
s = s + 2.0_real64 * f(a + real(i, real64) * h)
else
s = s + 4.0_real64 * f(a + real(i, real64) * h)
end if
end do
s = s * h / 3.0_real64
end function simpson
pure real(real64) function gauss2(f, a, b) result(s)
procedure(integrand) :: f
real(real64), intent(in) :: a, b
real(real64) :: mid, half
real(real64), parameter :: node = 0.5773502691896258_real64
mid = 0.5_real64 * (a + b)
half = 0.5_real64 * (b - a)
s = half * (f(mid - half * node) + f(mid + half * node))
end function gauss2
pure real(real64) function gauss5(f, a, b) result(s)
procedure(integrand) :: f
real(real64), intent(in) :: a, b
real(real64), parameter :: nodes(5) = [0.0_real64, &
0.5384693101056831_real64, -0.5384693101056831_real64, &
0.9061798459386640_real64, -0.9061798459386640_real64]
real(real64), parameter :: weights(5) = [0.5688888888888889_real64, &
0.4786286704993665_real64, 0.4786286704993665_real64, &
0.2369268850561891_real64, 0.2369268850561891_real64]
real(real64) :: mid, half
integer :: i
mid = 0.5_real64 * (a + b)
half = 0.5_real64 * (b - a)
s = 0.0_real64
do i = 1, 5
s = s + weights(i) * f(mid + half * nodes(i))
end do
s = s * half
end function gauss5
end module rules
The interface declares a pure function of one
real64 argument that returns a real64, and every
routine marks its function argument with procedure(integrand). The
pure attribute is not decoration: it lets the compiler inline the
call and forbids the function from modifying global state, which is exactly the
guarantee an integrator needs. The three rules differ only in the weights they
give the sample points.
program i1
use iso_fortran_env, only: real64
use rules
implicit none
integer, parameter :: steps(4) = [10, 100, 1000, 10000]
real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64)
integer :: i
print '(a)', ' n x^2 on [0,1] error sin on [0,pi] error'
do i = 1, size(steps)
print '(i5, f15.12, es12.3, f16.12, es12.3)', steps(i), &
trapezoid(square, 0.0_real64, 1.0_real64, steps(i)), &
trapezoid(square, 0.0_real64, 1.0_real64, steps(i)) - &
1.0_real64 / 3.0_real64, &
trapezoid(sine, 0.0_real64, pi, steps(i)), &
trapezoid(sine, 0.0_real64, pi, steps(i)) - 2.0_real64
end do
contains
pure real(real64) function square(x) result(y)
real(real64), intent(in) :: x
y = x * x
end function square
pure real(real64) function sine(x) result(y)
real(real64), intent(in) :: x
y = sin(x)
end function sine
end program i1
The integral of x**2 from zero to one is one third and the integral
of the sine over a half period is two, so every line shows the rule's error
next to its value. Each tenfold increase in n divides the error by
exactly one hundred, which is the signature of a second-order rule: the error
falls with the square of the step width. Ten thousand intervals still leave an
error of 1.6e-8 on the sine, so the trapezoid rule is adequate for a quick
estimate and wasteful for an accurate one.
$ gfortran rules.f90 i1.f90 -o i1 $ ./i1 n x^2 on [0,1] error sin on [0,pi] error 10 0.335000000000 1.667E-03 1.983523537509 -1.648E-02 100 0.333350000000 1.667E-05 1.999835503887 -1.645E-04 1000 0.333333500000 1.667E-07 1.999998355066 -1.645E-06 10000 0.333333335000 1.667E-09 1.999999983551 -1.645E-08
Midpoint and Simpson
The midpoint rule evaluates the function at the centre of each part instead of at its ends, and its error has the opposite sign but the same size, which makes the two rules a pair: their average is the trapezoid rule again. Simpson's rule takes the next step and fits a parabola through three consecutive points, giving the weights 1, 4, 1 to every pair of parts. The result is exact for any polynomial of degree three, and its error falls with the fourth power of the step, so a tenfold refinement divides it by ten thousand.
program i2
use iso_fortran_env, only: real64
use rules
implicit none
integer, parameter :: n = 8
real(real64) :: exact, e1, e2, m1, m2, s1, s2
exact = exp(1.0_real64) - 1.0_real64
e1 = trapezoid(exponential, 0.0_real64, 1.0_real64, n)
e2 = trapezoid(exponential, 0.0_real64, 1.0_real64, 2 * n)
m1 = midpoint(exponential, 0.0_real64, 1.0_real64, n)
m2 = midpoint(exponential, 0.0_real64, 1.0_real64, 2 * n)
s1 = simpson(exponential, 0.0_real64, 1.0_real64, n)
s2 = simpson(exponential, 0.0_real64, 1.0_real64, 2 * n)
print '(a, f12.10)', 'the exact integral: ', exact
print '(a)', 'rule n=8 error n=16 error ratio'
print '(a, 2es13.3, f9.2)', 'trapezoid', e1 - exact, e2 - exact, &
(e1 - exact) / (e2 - exact)
print '(a, 2es13.3, f9.2)', 'midpoint ', m1 - exact, m2 - exact, &
(m1 - exact) / (m2 - exact)
print '(a, 2es13.3, f9.2)', 'simpson ', s1 - exact, s2 - exact, &
(s1 - exact) / (s2 - exact)
print '(a, f16.12, 2x, f16.12)', 'simpson values: ', s1, s2
print '(a, f0.15)', 'simpson on x^3, 0 to 1: ', &
simpson(cube, 0.0_real64, 1.0_real64, 2)
print '(a, f0.15)', 'the exact value: ', 0.25_real64
contains
pure real(real64) function exponential(x) result(y)
real(real64), intent(in) :: x
y = exp(x)
end function exponential
pure real(real64) function cube(x) result(y)
real(real64), intent(in) :: x
y = x * x * x
end function cube
end program i2
With eight intervals the trapezoid and midpoint errors are 2.2e-3 and -1.1e-3,
and doubling the count divides both by four, exactly as the second-order theory
predicts. Simpson starts at 2.3e-6, three orders of magnitude smaller, and its
ratio is 15.98 rather than 16 because the next term of the expansion has begun
to contribute. The last two lines show the exactness property directly: two
intervals of Simpson's rule integrate x**3 to 0.25 without an error
in the last digit.
$ gfortran rules.f90 i2.f90 -o i2 $ ./i2 the exact integral: 1.7182818285 rule n=8 error n=16 error ratio trapezoid 2.237E-03 5.593E-04 4.00 midpoint -1.118E-03 -2.796E-04 4.00 simpson 2.326E-06 1.456E-07 15.98 simpson values: 1.718284154700 1.718281974052 simpson on x^3, 0 to 1: .250000000000000 the exact value: .250000000000000
Refining until the answer settles
In practice the exact value is unknown, so the usual stopping rule compares two
successive approximations: compute the integral with n parts,
double n, and stop when the change is below a tolerance. The change
is not the error, but for a converging rule it bounds the error to within a
small factor, and it costs one extra evaluation of the whole sum rather than a
proof. The loop below prints every round so that the fall of the change is
visible.
program i3
use iso_fortran_env, only: real64
use rules
implicit none
real(real64) :: coarse, fine, exact
integer :: n, rounds
exact = exp(1.0_real64) - 1.0_real64
n = 2
coarse = simpson(exponential, 0.0_real64, 1.0_real64, n)
rounds = 0
do
fine = simpson(exponential, 0.0_real64, 1.0_real64, 2 * n)
rounds = rounds + 1
print '(a, i6, a, f16.13, a, es10.2)', 'n = ', 2 * n, ', value ', &
fine, ', change ', abs(fine - coarse)
if (abs(fine - coarse) < 1.0e-13_real64) exit
if (n > 1000000) exit
coarse = fine
n = 2 * n
end do
print '(a, i0)', 'rounds: ', rounds
print '(a, f18.15)', 'the result: ', fine
print '(a, f18.15)', 'the exact value: ', exact
print '(a, es10.2)', 'the true error: ', abs(fine - exact)
contains
pure real(real64) function exponential(x) result(y)
real(real64), intent(in) :: x
y = exp(x)
end function exponential
end program i3
Ten rounds take the interval count from four to 2048, and the change between
successive approximations falls by a factor of about sixteen each time, which is
the fourth-order convergence of Simpson's rule in action. The final value agrees
with e - 1 to 1.3e-15, close to the limit of double precision. The
loop also has a hard stop after a million intervals, because a tolerance that is
too small would otherwise run until the count overflows.
$ gfortran rules.f90 i3.f90 -o i3 $ ./i3 n = 4, value 1.7183188419217, change 5.42E-04 n = 8, value 1.7182841546999, change 3.47E-05 n = 16, value 1.7182819740519, change 2.18E-06 n = 32, value 1.7182818375618, change 1.36E-07 n = 64, value 1.7182818290280, change 8.53E-09 n = 128, value 1.7182818284946, change 5.33E-10 n = 256, value 1.7182818284613, change 3.33E-11 n = 512, value 1.7182818284592, change 2.08E-12 n = 1024, value 1.7182818284591, change 1.35E-13 n = 2048, value 1.7182818284590, change 4.44E-15 rounds: 10 the result: 1.718281828459046 the exact value: 1.718281828459045 the true error: 1.33E-15
Bisection
Root finding starts from a bracket: two arguments at which the function has opposite signs, so that by the intermediate value theorem a root lies between them. Bisection halves the bracket at every step, keeping the half in which the sign change survives. Each step gains exactly one binary digit, so the number of steps to reach a tolerance depends only on the width of the first bracket and not on the function, which makes bisection slow but completely predictable.
module roots
use iso_fortran_env, only: real64
implicit none
private
public :: bisection, newton, secant
abstract interface
pure real(real64) function scalar_function(x) result(y)
import :: real64
real(real64), intent(in) :: x
end function scalar_function
end interface
contains
real(real64) function bisection(f, a, b, tolerance, steps) result(root)
procedure(scalar_function) :: f
real(real64), intent(in) :: a, b, tolerance
integer, intent(out) :: steps
real(real64) :: left, right, middle
left = a
right = b
steps = 0
do while (abs(right - left) > tolerance)
middle = 0.5_real64 * (left + right)
if (f(left) * f(middle) <= 0.0_real64) then
right = middle
else
left = middle
end if
steps = steps + 1
end do
root = 0.5_real64 * (left + right)
end function bisection
real(real64) function newton(f, df, guess, tolerance, steps) result(root)
procedure(scalar_function) :: f, df
real(real64), intent(in) :: guess, tolerance
integer, intent(out) :: steps
real(real64) :: x
x = guess
steps = 0
do while (abs(f(x)) > tolerance .and. steps < 50)
x = x - f(x) / df(x)
steps = steps + 1
end do
root = x
end function newton
real(real64) function secant(f, x0, x1, tolerance, steps) result(root)
procedure(scalar_function) :: f
real(real64), intent(in) :: x0, x1, tolerance
integer, intent(out) :: steps
real(real64) :: previous, current, following
previous = x0
current = x1
steps = 0
do while (abs(current - previous) > tolerance .and. steps < 50)
following = current - f(current) * (current - previous) / &
(f(current) - f(previous))
previous = current
current = following
steps = steps + 1
end do
root = current
end function secant
end module roots
The module keeps the three methods together and shares one interface. Bisection
needs no derivative, Newton's method needs the derivative as a second procedure,
and the secant method needs only the function. All three return the root and
report the number of steps through an intent(out) argument, which
is the standard way a Fortran routine returns a second result.
program i4
use iso_fortran_env, only: real64
use roots
implicit none
integer, parameter :: count = 5
real(real64), parameter :: tolerances(count) = [1.0e-2_real64, &
1.0e-4_real64, 1.0e-6_real64, 1.0e-9_real64, 1.0e-12_real64]
real(real64) :: root
integer :: i, steps
print '(a)', 'tolerance steps root f(root)'
do i = 1, count
root = bisection(cubic, 1.0_real64, 2.0_real64, tolerances(i), steps)
print '(es10.1, i8, f18.12, es13.3)', tolerances(i), steps, root, &
cubic(root)
end do
print '(a, f0.12)', 'x^3 - x - 2 at the root: ', cubic(1.521379706804_real64)
contains
pure real(real64) function cubic(x) result(y)
real(real64), intent(in) :: x
y = x * x * x - x - 2.0_real64
end function cubic
end program i4
The equation x**3 - x - 2 = 0 has one real root, near 1.5214, and
the bracket from one to two contains it because the function is -2 at the left
end and 4 at the right. Tightening the tolerance by a factor of a thousand adds
about ten steps every time, the ten binary digits of that factor, and the final
root is correct to twelve decimals. The value of the function at that point is
-3e-12 rather than zero, which is as close as double precision allows for a
number of this size.
$ gfortran roots.f90 i4.f90 -o i4 $ ./i4 tolerance steps root f(root) 1.0E-02 7 1.519531250000 -1.097E-02 1.0E-04 14 1.521392822266 7.796E-05 1.0E-06 20 1.521379947662 1.432E-06 1.0E-09 30 1.521379706915 6.592E-10 1.0E-12 40 1.521379706805 2.359E-12 x^3 - x - 2 at the root: -.000000000003
Newton's method
Newton's method replaces the curve by its tangent at the current guess and jumps
to where the tangent crosses zero, x - f(x)/f'(x). Near a simple
root the error is squared at every step, so the number of correct digits doubles:
three steps take an error of a hundredth to one of ten to the minus six, and the
fourth reaches the limit of the type. The price is the derivative, which must be
written by hand or approximated, and a good starting point, since a tangent can
throw the iterate far away.
program i5
use iso_fortran_env, only: real64
use roots
implicit none
real(real64) :: x, following, root
integer :: k, steps
print '(a)', ' k x(k) step f(x)'
x = 1.5_real64
do k = 1, 6
following = x - cubic(x) / derivative(x)
print '(i3, f20.15, es12.3, es12.3)', k, following, &
abs(following - x), cubic(following)
x = following
end do
root = newton(cubic, derivative, 1.5_real64, 1.0e-14_real64, steps)
print '(a, i0)', 'newton needed ', steps
print '(a, f0.15)', 'the root: ', root
root = newton(cubic, derivative, 10.0_real64, 1.0e-14_real64, steps)
print '(a, i0, a, f0.15)', 'from 10 it needed ', steps, ', root ', root
contains
pure real(real64) function cubic(x) result(y)
real(real64), intent(in) :: x
y = x * x * x - x - 2.0_real64
end function cubic
pure real(real64) function derivative(x) result(y)
real(real64), intent(in) :: x
y = 3.0_real64 * x * x - 1.0_real64
end function derivative
end program i5
The table shows the step lengths 2.2e-2, 3.6e-4, 9.9e-8 and 7.6e-15: each one is roughly the square of the one before, the fingerprint of quadratic convergence. The routine needs four steps from 1.5, and the fifth and sixth change nothing because the root is already exact to every stored digit. Starting from 10 instead needs ten steps, which shows that the method converges quickly only near the root; a program that cannot supply a good guess combines bisection with Newton, using bisection to get close and Newton to finish.
$ gfortran roots.f90 i5.f90 -o i5 $ ./i5 k x(k) step f(x) 1 1.521739130434783 2.174E-02 2.137E-03 2 1.521379805964786 3.593E-04 5.894E-07 3 1.521379706804575 9.916E-08 4.530E-14 4 1.521379706804568 7.550E-15 0.000E+00 5 1.521379706804568 0.000E+00 0.000E+00 6 1.521379706804568 0.000E+00 0.000E+00 newton needed 4 the root: 1.521379706804568 from 10 it needed 10, root 1.521379706804568
The secant method
The secant method removes the derivative by approximating it with the slope
between the last two iterates, (f(x1) - f(x0)) / (x1 - x0), which
is the secant line through the two points. It therefore needs two starting
values instead of one and no derivative at all, and its order of convergence is
the golden ratio, about 1.618, slower than Newton's square but faster than
bisection's halving.
program i6
use iso_fortran_env, only: real64
use roots
implicit none
real(real64) :: x0, x1, following, root
integer :: k, steps
print '(a)', ' k x(k) step f(x)'
x0 = 1.0_real64
x1 = 2.0_real64
do k = 1, 7
following = x1 - cubic(x1) * (x1 - x0) / (cubic(x1) - cubic(x0))
print '(i3, f20.15, es12.3, es12.3)', k, following, &
abs(following - x1), cubic(following)
x0 = x1
x1 = following
end do
root = secant(cubic, 1.0_real64, 2.0_real64, 1.0e-14_real64, steps)
print '(a, i0, a, f0.15)', 'secant needed ', steps, ', root ', root
root = newton(cubic, derivative, 1.0_real64, 1.0e-14_real64, steps)
print '(a, i0, a, f0.15)', 'newton needed ', steps, ', root ', root
contains
pure real(real64) function cubic(x) result(y)
real(real64), intent(in) :: x
y = x * x * x - x - 2.0_real64
end function cubic
pure real(real64) function derivative(x) result(y)
real(real64), intent(in) :: x
y = 3.0_real64 * x * x - 1.0_real64
end function derivative
end program i6
The step lengths fall as 0.67, 0.13, 0.068, 0.010, 4.5e-4, 3.4e-6 and 1.2e-9, which is neither halving nor squaring but the golden ratio at work. To reach the same tolerance the secant method needs eight steps against Newton's six, and it evaluates the function once per step instead of twice, since the old value is reused. When the derivative is expensive or awkward, the secant method is usually the better trade.
$ gfortran roots.f90 i6.f90 -o i6 $ ./i6 k x(k) step f(x) 1 1.333333333333333 6.667E-01 -9.630E-01 2 1.462686567164179 1.294E-01 -3.333E-01 3 1.531169432141204 6.848E-02 5.863E-02 4 1.520926420515280 1.024E-02 -2.693E-03 5 1.521376316669744 4.499E-04 -2.015E-05 6 1.521379707984872 3.391E-06 7.015E-09 7 1.521379706804564 1.180E-09 -1.799E-14 secant needed 8, root 1.521379706804568 newton needed 6, root 1.521379706804568
Gaussian quadrature
All the rules so far sample the function at equally spaced points. Gaussian
quadrature chooses the points and the weights together so that a rule with
n points integrates every polynomial of degree up to
2n - 1 exactly, which is the best any n-point rule can
do. For two points the nodes are at plus and minus one over the square root of
three with equal weights, and for five points the nodes and weights are the
tabulated constants below. The interval is mapped to the standard one from -1
to 1 by a linear change of variable.
program i7
use iso_fortran_env, only: real64
use rules
implicit none
real(real64) :: exact
exact = exp(1.0_real64) - 1.0_real64
print '(a, f16.13)', 'trapezoid, n = 100: ', &
trapezoid(exponential, 0.0_real64, 1.0_real64, 100)
print '(a, f16.13)', 'simpson, n = 100: ', &
simpson(exponential, 0.0_real64, 1.0_real64, 100)
print '(a, f16.13)', 'gauss, two points: ', gauss2(exponential, &
0.0_real64, 1.0_real64)
print '(a, f16.13)', 'gauss, five points: ', gauss5(exponential, &
0.0_real64, 1.0_real64)
print '(a, f16.13)', 'the exact value: ', exact
print '(a, es10.2)', 'trapezoid error: ', &
abs(trapezoid(exponential, 0.0_real64, 1.0_real64, 100) - exact)
print '(a, es10.2)', 'simpson error: ', &
abs(simpson(exponential, 0.0_real64, 1.0_real64, 100) - exact)
print '(a, es10.2)', 'gauss two points: ', &
abs(gauss2(exponential, 0.0_real64, 1.0_real64) - exact)
print '(a, es10.2)', 'gauss five points: ', &
abs(gauss5(exponential, 0.0_real64, 1.0_real64) - exact)
print '(a, f0.15)', 'gauss5 on x^7: ', gauss5(seventh, &
0.0_real64, 1.0_real64)
print '(a, f0.15)', 'the exact 1/8: ', 0.125_real64
print '(a, f0.15)', 'simpson on x^7, n=2:', simpson(seventh, &
0.0_real64, 1.0_real64, 2)
contains
pure real(real64) function exponential(x) result(y)
real(real64), intent(in) :: x
y = exp(x)
end function exponential
pure real(real64) function seventh(x) result(y)
real(real64), intent(in) :: x
y = x**7
end function seventh
end program i7
Five points reach an error of 6.5e-13 on the exponential, while Simpson's rule
with a hundred intervals, that is two hundred and one evaluations, is stuck at
9.6e-11. The two-point rule is worse than Simpson here because the exponential
is not a polynomial, but it costs two evaluations. The last three lines show the
exactness property: five points integrate x**7 to 0.125, while two
intervals of Simpson's rule give 0.171875 for the same integral.
$ gfortran rules.f90 i7.f90 -o i7 $ ./i7 trapezoid, n = 100: 1.7182961474504 simpson, n = 100: 1.7182818285545 gauss, two points: 1.7178963780075 gauss, five points: 1.7182818284584 the exact value: 1.7182818284590 trapezoid error: 1.43E-05 simpson error: 9.55E-11 gauss two points: 3.85E-04 gauss five points: 6.53E-13 gauss5 on x^7: .125000000000000 the exact 1/8: .125000000000000 simpson on x^7, n=2:.171875000000000
Next steps
The procedure arguments article
explains the abstract interface and the procedure dummy arguments
that make these routines generic, and the
floating-point article explains the
rounding that limits the final digits and the cancellation that a difference
such as f(x1) - f(x0) can suffer. The
intrinsic numeric functions article
covers abs and exp, the
arrays article covers the array constants of the
five-point rule, and the loops article explains
the do while construct that drives every iteration.
Source
Fortran 2018 standard (J3/18-007r1)
In this article we have implemented the trapezoid, midpoint and Simpson rules with a function passed as a procedure argument, measured their second- and fourth-order convergence on integrals with known values, refined an integral until successive answers agreed, found a cubic's root by bisection and watched the step count grow with the tolerance, used Newton's method and saw quadratic convergence break down far from the root, solved the same equation by the secant method without a derivative, and integrated an exponential and a seventh power with two- and five-point Gaussian quadrature. The previous article, Fortran error handling, supplies the status codes that a numerical routine returns when it cannot converge.
Author
List all Fortran tutorials.