Numerical Integration
Ward Cheney and David Kincaid: Chapter 5
Computer Problems 5-1: 1,2,3,5,6
Computer Problems 5-2: 1 →4, 8,10
Numerical integration Ch5 1
Definite and Indefinite Integrals
function
Number
Lower and Upper Sums: The existence of the definite integral of a nonnegative
function f on a closed interval [a, b] is based on an interpretation of that integral as
the area under the graph of f. The definite integral is defined by means of two
concepts, the lower sums of f and the upper sums of f ; these are approximations
to the area under the curve.
Let P be a partition of the interval [a , b ] given by
with partition points x0, x1, x2,..., xn that divide the interval [a,b] into n subintervals
[xi, xi+1]. Now denote by mi the lower bound (infimum or inf) of f (x) on the
subinterval [xi, xi+1]. In symbols,
Numerical integration Ch5 2
Lower and Upper Sums
Similarly, the upper bound (supremum or sup) of f (x) on the subinterval
[xi, xi+1] is denoted by Mi
The lower sums and upper sums of f corresponding to the given
partition P are defined as
Numerical integration Ch5 3
Lower and Upper Sums
f(x) is increasing on [0, 1], the least (greatest) value of f on the subinterval [xi, xi+1] ] occurs at xi (xi+1)
It is intuitively clear that the upper sum overestimates the area under the curve, and
the lower sum underestimates it. OneNumerical
canintegration
showCh5Error 4
sums-uplow.f95
program sums
integer :: i
For partitions, we take equally spaced points in [0 , 1]. Thus,
real :: h, sum, suml, sumu, x if there are to be n subintervals
integer, parameter :: n = 1000
real, parameter :: a = 0.0 ; b = 1.0
here
h = 1.0/real(n)
sum = 0.0 f(x) is decreasing on [0, 1], the least value of f on the
do i = 1,n
x = a + real(i)*h
subinterval [xi, xi+1] ] occurs at xi+1 . Similarly, the greatest
sum = sum + f(x) value occurs at xi. Hence, mi = f (xi+1) and Mi = f (xi)
end do
suml = h*sum
sumu = suml + h*(f(a) - f(b))
print *,"suml =",suml,"sumu =",sumu
end program sums
function f(x)
real, intent(in) :: x As the formulas are almost the same, we can write for
f = 1.0/exp(x*x) computational benefits ( Homework assignment)
end function f
Numerical integration Ch5 5
Example
If the integral is to be computed with absolute error < 0.5×10−3 and if
we are going to use upper and lower sums with a uniform partition, how many
subintervals are needed?
The integrand, f (x)=exp(cos(x)), is a decreasing function on
the given interval [0,π] (Why ?) Therefore, we have
error of at most To meet the error criterion we must have
Note: all the middle terms
cancel out except o and pi.
Numerical integration Ch5 6
Trapezoid Rule
• Based on an estimation of the area beneath a curve using trapezoids.
• The estimation of is approached by first dividing the interval [a,b] into
subintervals according to the partition
• For each such partition of the interval an estimation of the integral by the
trapezoid rule is obtained.
Numerical integration Ch5 7
Trapezoid Rule
The area is equal to the base times the average height, and we have the basic
trapezoid rule for the subinterval [xi, xi+1] .
Therefore, total area = sum of all these areas, i.e.
This formula is called the composite trapezoid rule.
Numerical integration Ch5 8
Uniform Spacing (trapezoid.f95)
• In a uniform partition of the interval, the division points xi are equally spaced,
Think of h as the step size.
• Therefore, we have
• Computationally cheaper formula is (show it)
• Solve the integral numerically
Numerical integration Ch5 9
trapezoid.f95
1. Run the code trapezoid.f95 and observe the output and compare your
computed value with exact answer.
2. Change the size of step (h) or, equally to say that, change number of steps (n)
and see what happened to your predicted answers. Obviously, we expect that
as h →0, the relative error must also →0. (Check this)
3. Using the same program evaluate the integral
4. By changing different values of h show that in the trapezoidal rule the error
term is ~ O (h2)
5. The Error formula for Trapezoidal rule is (we don’t need to prove it)
Numerical integration Ch5 10
Example
How many points should be used while computing the following integral with an
error of at most 0.5× 10−4?
Homework: How many subintervals are needed to approximate
Numerical integration Ch5 11
Recursive Trapezoid Formula for Equal Subintervals
Composite trapezoid rule when the interval [a, b] is subdivided into 2n equal parts
Romberg algorithm: It denotes the result of applying composite trapezoid rule with
2n equal subintervals.
• In the Romberg algorithm, it will also be necessary to have a means of
computing R (n , 0) from R (n − 1, 0) w/o involving un-needed evaluations of f. If
R (n − 1, 0) has been computed and R (n , 0) is to be computed, we use the
identity
Numerical integration Ch5 (proof On next slide)
12
It is desirable to compute the bracketed expression with as little additional work as
possible. Fixing h = (b − a )/ 2n for the analysis and putting
r(0,0) in code
Then
Here, we have taken account of the fact that each term in the first sum that
corresponds to an even value of i is canceled by a term in the second sum. This
leaves only terms that correspond to odd values of i.
Numerical integration Ch5 13
Romberg method
You don’t need to derive all the following, but just use them to build your own
FORTRAN code to evaluate different definite Integrals.
The Romberg algorithm produces a triangular array of
numbers, all of which are numerical estimates of the
definite integral
If you still want to derive it see Section 4.3 of your book.
Numerical integration Ch5 14
Algorithm
Romberg integration is technique that can also be used not only for Trapezoidal method, but
also for Simpson and other integration methods, but here we only discussed it for Trapezoidal. It
is an extremely efficient way of improving the algorithms, and reducing the errors in different
integration methods, i.e., very high accuracy by using just few functions.
Numerical integration Ch5 15
romberg.f95
Romberg array for three separate functions:
1. f(x)= 1/(1+x)
2. g(x)=exp(x)
3. p(x)= sqrt(x)
subroutine rombrg(f,a,b,n,r)
m=1
real, dimension (0:n, 0:n) :: r
r(i,0) = 0.5*r(i-1,0) + h*sum
integer, intent (in) :: n
print "(a, i1, a, 5e22.14)", "r(", i,",0) =", r(i,0)
real, intent(in) :: a, b
do j = 1,i
real :: h, sum
m = 4*m
integer :: i, k, j
r(i,j) = r(i,j-1) + (r(i,j-1) - r(i-1,j-1))/real(m - 1)
h=b-a
print "(a, i1, a, i1, a, 5e22.14)", "r(", i, ",", j,") = ", r(i,j)
r(0,0) = 0.5*h*(f(a) + f(b))
end do
do i = 1, n
end do
h = 0.5*h
end subroutine rombrg
sum = 0.0
do k = 1,2**i -1,2
sum = sum + f(a + h*real(k))
end do
Numerical integration Ch5 16