# 2.12 NUMERICAL INTEGRATION

```#
#
# WORKSHEET # 23
#
#  NUMERICAL INTEGRATION
#
#
--------------------------------------------------------------------------------
# In this worksheet we look at several rules for numerical evaluation of integrals
# Int(f(x),x=a..b).  The interval [a,b] is subdivided into n equal sub-intervals, all of length
# (b-a)/n.  The end points of the subintervals are denoted by X(i,n).
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
> h:=n->(b-a)/n;

b - a
h := n -> -----
n
--------------------------------------------------------------------------------
> X:=(i,n)->a+i*(b-a)/n;

i (b - a)
X := (i,n) -> a + ---------
n
--------------------------------------------------------------------------------
# The mid point rule is given as follows.  The interval [a,b] is subdivided into n equal sub-
# intervals, all of length h(n), and then the function f(x) is evaluated at each of the mid
# points (1/2)*(X(i,n)+X(i+1,n)), i=0..n-1.
--------------------------------------------------------------------------------
> Mid:=n->h(n)*sum(f((X(i,n)+X(i+1,n))/2),i=0..n-1);

/n - 1                                 \
|-----                                 |
| \                                    |
Mid := n -> h(n) |  )   f(1/2 X(i, n) + 1/2 X(i + 1, n))|
| /                                    |
|-----                                 |
\i = 0                                 /
--------------------------------------------------------------------------------
# Notice that there are n functional evaluations in the mid point rule and that the product
# of the multiplier h(n)=(b-a)/n and the number of functional evaluations is b-a.
--------------------------------------------------------------------------------
# In the trapezoidal rule we also subdivide the interval [a,b] into n equal parts, but we
# approximate the integral in each subinterval by a trapezoid.
--------------------------------------------------------------------------------
> Trap:=n->(h(n)/2)*sum(f(X(i,n))+f(X(i+1,n)),i=0..n-1);

/n - 1                              \
|-----                              |
| \                                 |
Trap := n -> 1/2 h(n) |  )   (f(X(i, n)) + f(X(i + 1, n)))|
| /                                 |
|-----                              |
\i = 0                              /
--------------------------------------------------------------------------------
# Again the product of the multiplier h(n)/2 and the number of functional evaluations is
# b-a.  Let M[n] and T[n] denote the mid point and trapezoidal rules.  If f(x) is twice
# continuously differentiable on the interval [a,b] then the error formulas are given as
# follows.  The second order derivatives are calculated at some unknown points t, u in the
# interval [a,b].
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=M[n]+D(D(f))(t)*(b-a)^3/(24*n^2);

b
/                        (2)              3
|                        D   (f)(t) (b - a)
|  f(x) dx = M[n] + 1/24 -------------------
|                                  2
/                                  n
a
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[n]-D(D(f))(u)*(b-a)^3/(12*n^2);

b
/                        (2)              3
|                        D   (f)(u) (b - a)
|  f(x) dx = T[n] - 1/12 -------------------
|                                  2
/                                  n
a
--------------------------------------------------------------------------------
# Notice that if we ignore the second order derivatives in the error formulas, or if t is near
# u,  then the errors in these formulas have opposite signs and the mid point rule is twice as
# accurate as the trapezoidal rule.
--------------------------------------------------------------------------------
# In order to get an idea how accurate these rules are we apply them to some integral where
# we already know the answer, namely Int(1/(x+1),x=0..1).
--------------------------------------------------------------------------------
> a:=0:b:=1:
--------------------------------------------------------------------------------
> f:=x->1/(x+1);

1
f := x -> -----
x + 1
--------------------------------------------------------------------------------
> J:=int(f(x),x=a..b);

J := ln(2)
--------------------------------------------------------------------------------
> evalf(ln(2));

.6931471806
--------------------------------------------------------------------------------
> seq([n,evalf(J-Mid(n)),evalf(J-Trap(n))],n=1..5);

[1, .0264805139, -.0568528194], [2, .0074328949, -.0151861527],

[3, .0033924908, -.0068528194], [4, .0019272894, -.0038766289],

[5, .0012392949, -.0024877400]
--------------------------------------------------------------------------------
# We expect the error J-Mid(n) to be positive since, according to the above error formula,
# J-Mid(n)=2*(t+1)^{-2}/24*n^2 for some t between 0 and 1.  Likewise we expect the
# error for the trapezoidal rule to be negative.  The above calculations corroborate this.
--------------------------------------------------------------------------------
# The next calculation reveals that as n gets larger the errors in the mid point and
# trapezoidal rules are behaving like constants times n^2.  from the graph it would seem
# that Limit(n^2*(J-Mid(n),n=infinity) is approximately 0.03, whereas
# Limit(n^2*(J-Trap(n),n=infinity) is approximately -0.06.  Notice the difference in signs
# and the fact that one is double the other in absolute value.
--------------------------------------------------------------------------------
> plot({[seq([n,n^2*evalf(J-Mid(n))],n=1..20)],[seq([n,n^2*evalf(J-Trap(n))],n=1..20)]},style=point);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# If the function f(x) is not continuously differentiable on the interval then the error
# formulas may no longer be valid.  Consider the next example, Int(sqrt(x),x=0..1).
--------------------------------------------------------------------------------
> f:=x->sqrt(x);

f := sqrt
--------------------------------------------------------------------------------
> J:=int(f(x),x=a..b);

J := 2/3
--------------------------------------------------------------------------------
> plot({[seq([n,n^2*evalf(J-Mid(n))],n=1..20)],[seq([n,n^2*evalf(J-Trap(n))],n=1..20)]},style=point);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# These graphs are not converging to a horizontal asymptote.
--------------------------------------------------------------------------------
> f:='f':
--------------------------------------------------------------------------------
# The error formulas for the mid point and trapezoidal rules suggest that the linear
# combination (2/3)M[n]+(1/3)T[n] is a better approximation, because if t, u are reasonably
# close then this combination should reduce the error.  Note that if t=u then this
# combination completely eliminates the error.  Simpson's rule is just this combination of
# the mid point and trapezoidal rules.  After some algebra and rewriting the rule becomes:
--------------------------------------------------------------------------------
> Simpson:=n->(h(n)/3)*sum(f(X(2*i,n))+4*f(X(2*i+1,n))+f(X(2*i+2,n)),i=0..n/2-1);

Simpson :=

/1/2 n - 1                                                       \
|  -----                                                         |
|   \                                                            |
n -> 1/3 h(n) |    )     (f(X(2 i, n)) + 4 f(X(2 i + 1, n)) + f(X(2 i + 2, n)))|
|   /                                                            |
|  -----                                                         |
\  i = 0                                                         /
--------------------------------------------------------------------------------
# In this case the interval [a,b] is subdivided into n equal parts, where n must be even.  The
# number of functional evaluations is 3n and again the product with the multiplier h(n)/3 is
# b-a, as it should be. Set S[n] equal to Simpsons rule.  If f(x) is 4 times continuously
# differentiable on the interval [a,b] then the error in Simpson's rule is given as follows,
# where t is just some unknown point in [a,b].
--------------------------------------------------------------------------------
> a:='a':b:='b':n:='n':f:='f':Int(f(x),x=a..b)=S[n]-(D@@4)(f)(t)*(b-a)^5/(180*n^4);

b
/                         (4)              5
|                         D   (f)(t) (b - a)
|  f(x) dx = S[n] - 1/180 -------------------
|                                   4
/                                   n
a
--------------------------------------------------------------------------------
# Let's apply Simpson's rule to Int(1/(x+1),x=0..1) in order to get an idea of how accurate
# this rule is.
--------------------------------------------------------------------------------
> a:=0:b:=1:
--------------------------------------------------------------------------------
> f:=x->1/(x+1);

1
f := x -> -----
x + 1
--------------------------------------------------------------------------------
> J:=ln(2);

J := ln(2)
--------------------------------------------------------------------------------
> seq([2*n,evalf(J-Simpson(2*n))],n=1..5);

-5
[2, -.0012972638], [4, -.0001067877], [6, -.0000226126], [8, -.73501*10  ],

-5
[10, -.30501*10  ]
--------------------------------------------------------------------------------
# We expect the error in Simpson's rule to be negative since the fourth order derivative of
# 1/(x+1) is 24*(x+1)^{-5}, which is positive on the interval [0,1].  In fact the error is given
# by  J-Simpson(2*n)=-24*(t+1)^{-5}/180*(2*n)^4, some t between 0 and 1.
--------------------------------------------------------------------------------
# We also expect (2*n)^{4}*( J-Simpson(2*n)) to approach a negative constant as n tends
# to infinity.
--------------------------------------------------------------------------------
> plot([seq([2*n,(2*n)^4*evalf(J-Simpson(2*n))],n=1..20)],style=point);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# An immediate consequence of the error formulas for the mid point and trapezoidal rules is
# that they give exact answers for polynomials of degree less than or equal to 1.  On the
# other hand, Simpson's rule gives the exact value for a polynomial whose degree is no
# more than 3.  There are other rules that are exact on polynomials of higher degree.
--------------------------------------------------------------------------------
> poly:=n->unapply(sum(c[j]*x^j,j=0..n),x);

n
-----
\          j
poly := n -> unapply(  )   c[j] x , x)
/
-----
j = 0
--------------------------------------------------------------------------------
> f:=poly(5);

2         3         4         5
f := x -> c[0] + c[1] x + c[2] x  + c[3] x  + c[4] x  + c[5] x
--------------------------------------------------------------------------------
> J:=f->int(f(x),x=a..b);

b
/
|
J := f ->  |  f(x) dx
|
/
a
--------------------------------------------------------------------------------
> J(f)-Mid(1);

11          13
1/12 c[2] + 1/8 c[3] + ---- c[4] + ---- c[5]
80          96
--------------------------------------------------------------------------------
> J(f)-Trap(1);

- 1/6 c[2] - 1/4 c[3] - 3/10 c[4] - 1/3 c[5]
--------------------------------------------------------------------------------
> J(f)-Simpson(2);

- 1/120 c[4] - 1/48 c[5]
--------------------------------------------------------------------------------
# The end points a, b do not appear in the answers since we previously assigned the values
# 0, 1 to them.   We do not see the coefficients c[0] or c[1] appearing in the first 2 answers
# because the mid point and trapezoidal rules are exact on polynomials of degree 1, and we
# do not see c[0], c[1], c[2] or c[3] in J(f)-Simpson(2) because Simpson's rule is exact on
# polynomials of degree 3.  In the next few calculations we rederive Simpson's rule by
# searching for a rule which is exact on polynomials of degree less than or equal to 3.  This
# is just a warm-up exercise  for what is to come.
--------------------------------------------------------------------------------
> rule:=n->h(n)*sum(d[0]*f(X(2*i,n))+d[1]*f(X(2*i+1,n))+d[2]*f(X(2*i+2,n)),i=0..n/2-1);

rule := n -> h(n)

/1/2 n - 1                                                                    \
|  -----                                                                      |
|   \                                                                         |
|    )     (d[0] f(X(2 i, n)) + d[1] f(X(2 i + 1, n)) + d[2] f(X(2 i + 2, n)))|
|   /                                                                         |
|  -----                                                                      |
\  i = 0                                                                      /
--------------------------------------------------------------------------------
# Here we subdivide the interval [a,b] into n equal parts, where n is even, and then take a
# certain linear combination of functional values.  We want (d[0]+d[2]+d[2])*(n/2)*h(n) to
# be equal to b-a, and we want to determine the unknowns d[0], d[1] and d[2] so that this
# rule is exact on polynomials of degree 3.  It is sufficient to determine the d[j] so that
# rule(2) is exact.
#
# EXERCISE:  Why is this true?
--------------------------------------------------------------------------------
> J(f)-rule(2);

c[0] + 1/2 c[1] + 1/3 c[2] + 1/4 c[3] + 1/5 c[4] + 1/6 c[5] - 1/2 d[0] c[0]

- 1/2 d[1] (c[0] + 1/2 c[1] + 1/4 c[2] + 1/8 c[3] + 1/16 c[4] + 1/32 c[5])

- 1/2 d[2] (c[0] + c[1] + c[2] + c[3] + c[4] + c[5])
--------------------------------------------------------------------------------
# Can we choose the d[j]'s so that the coefficients of c[0], c[1], c[2] and c[3] are all zero?  In
# the next calculation we collect these coefficients and make them into equations.
--------------------------------------------------------------------------------
> eqns:={seq(coeff(expand("),c[j],1),j=0..3)};

eqns := {1/3 - 1/8 d[1] - 1/2 d[2], 1/2 - 1/4 d[1] - 1/2 d[2],

1/4 - 1/16 d[1] - 1/2 d[2], 1 - 1/2 d[0] - 1/2 d[1] - 1/2 d[2]}
--------------------------------------------------------------------------------
> solve(eqns);

{d[0] = 1/3, d[2] = 1/3, d[1] = 4/3}
--------------------------------------------------------------------------------
# These values of d[j] are exactly the values for Simpson's rule.
--------------------------------------------------------------------------------
> f:='f':
--------------------------------------------------------------------------------
# We can generalize the idea in this last calculation to come up with rules that are exact for
# all polynomials up to a certain degree.  These are called Newton-Cotes rules.  Such rules
# become more accurate as the degree of the polynomial goes up.  For example, one
# application of the fourth order Newton-Cotes rule to Int(f(x),x=a..b) is given by the
# following formula, where h=(b-a)/4 and f[j]=a+j*h, j=0..4.  We derive this formula below.
--------------------------------------------------------------------------------
> a:='a':b:='b':
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=(2*h/45)*(7*f[0]+32*f[1]+12*f[2]+32*f[3]+7*f[4])-(8*h^7/945)*(D@@6)(f)(t);

b
/
|
|  f(x) dx =
|
/
a

7  (6)
2/45 h (7 f[0] + 32 f[1] + 12 f[2] + 32 f[3] + 7 f[4]) - 8/945 h  D   (f)(t)
--------------------------------------------------------------------------------
> a:=0;b:=1;

a := 0

b := 1
--------------------------------------------------------------------------------
# To find higher order Newton-Cotes rules we choose positive integers n, k so that k divides
# n.  Then we assume that the following rule is exact on polynomials of high degree (degree
# k if k is odd and degree k+1 if k is even) and solve for the coefficients d[j].
--------------------------------------------------------------------------------
> rule:=(n,k)->h(n)*sum(sum(d[j]*f(X(k*i+j,n)),j=0..k),i=0..n/k-1);

/n/k - 1 /  k                        \\
| -----  |-----                      ||
|  \     | \                         ||
rule := (n,k) -> h(n) |   )    |  )   d[j] f(X(k i + j, n))||
|  /     | /                         ||
| -----  |-----                      ||
\ i = 0  \j = 0                      //
--------------------------------------------------------------------------------
# The interval [a,b] has been subdivided into n equal parts all of length h(n)=(b-a)/n.
# There are really n/k applications of the rule, one for each group of k+1 consecutive points
# X(k*i,n), X(k*i+1,n),...,X(k*i+k,n), i=0..n/k-1.  Notice that the coefficients for each group
# of points repeats, i.e., d[0], d[1],...,d[k], and recall that we have previously assigned a=0,
# b=1.
--------------------------------------------------------------------------------
> i:='i':j:='j':
--------------------------------------------------------------------------------
> rule(2,2);

1/2 d[0] f(0) + 1/2 d[1] f(1/2) + 1/2 d[2] f(1)
--------------------------------------------------------------------------------
> rule(3,3);

1/3 d[0] f(0) + 1/3 d[1] f(1/3) + 1/3 d[2] f(2/3) + 1/3 d[3] f(1)
--------------------------------------------------------------------------------
> rule(4,4);

1/4 d[0] f(0) + 1/4 d[1] f(1/4) + 1/4 d[2] f(1/2) + 1/4 d[3] f(3/4)

+ 1/4 d[4] f(1)
--------------------------------------------------------------------------------
> rule(k,k);

k
-----
\
)   d[j] f(j/k)
/
-----
j = 0
-----------------
k
--------------------------------------------------------------------------------
> rule(4,2);

1/4 d[0] f(0) + 1/4 d[1] f(1/4) + 1/4 d[2] f(1/2) + 1/4 d[0] f(1/2)

+ 1/4 d[1] f(3/4) + 1/4 d[2] f(1)
--------------------------------------------------------------------------------
> eqns:=k->eval({seq(subs(f=unapply(x^i,x),rule(k,k)=J(f)),i=0..k)});
Warning, `i` is implicitly declared local

eqns := k -> local i;

i
eval({seq(subs(f = unapply(x , x), rule(k, k) = J(f)), i = 0 .. k)})
--------------------------------------------------------------------------------
> eqns(2);

{1/8 d[1] + 1/2 d[2] = 1/3, 1/4 d[1] + 1/2 d[2] = 1/2,

1/2 d[0] + 1/2 d[1] + 1/2 d[2] = 1}
--------------------------------------------------------------------------------
> eqns(4);

{1/16 d[1] + 1/8 d[2] + 3/16 d[3] + 1/4 d[4] = 1/2,

1/4 d[0] + 1/4 d[1] + 1/4 d[2] + 1/4 d[3] + 1/4 d[4] = 1,

1/64 d[1] + 1/16 d[2] + 9/64 d[3] + 1/4 d[4] = 1/3,

27
1/256 d[1] + 1/32 d[2] + --- d[3] + 1/4 d[4] = 1/4,
256

81
1/1024 d[1] + 1/64 d[2] + ---- d[3] + 1/4 d[4] = 1/5}
1024
--------------------------------------------------------------------------------
> solve(");

14           64           14           64
{d[4] = ----, d[3] = ----, d[0] = ----, d[1] = ----, d[2] = 8/15}
45           45           45           45
--------------------------------------------------------------------------------
> subs(",rule(3,3));

14         64                         64
--- f(0) + --- f(1/3) + 8/45 f(2/3) + --- f(1)
135        135                        135
--------------------------------------------------------------------------------
> eval(subs(f=poly(4),J(f)-"));

31         689        203         331         4291
- --- c[0] - ---- c[4] - --- c[1] - ---- c[2] - ----- c[3]
135        2187        810        1215        14580
--------------------------------------------------------------------------------
> subs(solve(eqns(4)),rule(4,4));

16                          16
7/90 f(0) + ---- f(1/4) + 2/15 f(1/2) + ---- f(3/4) + 7/90 f(1)
45                          45
--------------------------------------------------------------------------------
# This last formula is Newton-Cotes rule of order 4, but it is actually exact for polynomials
# of degree 5.
--------------------------------------------------------------------------------

```