#
#
# 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.
--------------------------------------------------------------------------------
