Introduction

Topics in this Section

▸ Approximation of Derivatives if function values are knows.
▸ Higher order approximation.
▸ Approximation of Second order derivatives.
▸ Numerical Differentiation using Lagrange Interpolation
▸ Roundoff Errors and Numerical Differentialtions

2 / 33
Recall

Theorem (Taylor’s Theorem)

Assume that f ∈ C n [a, b], f (n+1) exists on [a, b] , and x0 ∈ [a, b].
Then for every x ∈ [a, b], there exists a number x1 (x) such that,

f (n)
f (x) = f (x0 ) + f ′ (x0 )(x − x0 ) + ⋯ + (x − x0 )n + Rn
n!
f (n+1) (x1 )
where Rn = (n+1)! (x − x0 )n+1 .

3 / 33
First order Derivative
▸ Derivative of a function:
f (x0 + h) − f (x0 )
f ′ (x0 ) = lim
h→0 h
if the limit exists, otherwise not defined.
▸ Assume f to be C 2 (a, b),then using Taylor’s Theorem:

h2 ′′
f (x0 + h) = f (x0 ) + hf ′ (x0 ) + f (ξ)
2
for some ξ.
▸ Rearranging terms:

f (x0 + h) − f (x0 ) h ′′
f ′ (x0 ) = − f (ξ).
h 2

4 / 33
First order Derivative: Two Point Formulas

▸ Forward Difference:
f (x0 + h) − f (x0 )
f ′ (x0 ) ≈ ,
h
with error
∣h∣ ′′
∣f (ξ)∣.
2
▸ Similarly, we have Backward Difference:

f (x0 ) − f (x0 − h)
f ′ (x0 ) ≈ ,
h
with error
∣h∣ ′′
∣f (ξ)∣.
2

5 / 33
First order Derivative: Three Point Formulas

▸ Assume f to be C 3 (a, b),then using Taylor’s Theorem:

h2 ′′ h3
f (x0 + h) = f (x0 ) + hf ′ (x0 ) + f (x0 ) + f (3) (ξ1 )
2 3!
for some ξ1 .
▸ Similarly,

h2 ′′ h3
f (x0 − h) = f (x0 ) − hf ′ (x0 ) + f (x0 ) − f (3) (ξ2 )
2 3!
for some ξ2 .

6 / 33
First order Derivative: Three Point Formulas

▸ Solving Both for f ′ (x0 ),we get,

f (x0 + h) − f (x0 − h) h2 (3)

f ′ (x0 ) = − f (ξ)
2h 6
where ξ is such that,

h2 (3) 1 h3 h3
f (ξ) = ( f (3) (ξ2 ) + f (3) (ξ1 ))
6 2 3! 3!

by Intermediate Value Theorem.

7 / 33
First order Derivative: Three Point Formulas

▸ Central Difference Formula:

f (x0 + h) − f (x0 − h)
f ′ (x0 ) ≈ ,
2h
with error
∣h∣2 3
∣f (ξ)∣.
6
for some ξ in (x0 − h, x0 + h)

8 / 33
First order Derivative: Three Point Formulas

▸ Similarly we can also get One sided Three point formula:

1
f ′ (x0 ) ≈ (−3f (x0 ) + 4f (x0 + h) − f (x0 + 2h)) ,
2h
with error
h2 (3)
∣f (ξ)∣
3
for some ξ in (x0 , x0 + 2h).
▸ Both three point formulas are second order accurate.
▸ Replacing h with −h we can get other side formula.

9 / 33
First order Derivative: Five-Point Formulas

▸ Central formula:
1
f ′ (x0 ) ≈ (f (x0 − 2h) − 8f (x0 − h) + 8f (x0 + h) − f (x0 + 2h)) ,
12h
with error
h4 (5)
∣f (ξ)∣
30
for some ξ in (x0 − 2h, x0 + 2h).
▸ Exercise: One sided five point formula.

10 / 33
Second order Derivative: Three Point Formulas

▸ Note that,

h4 (4)
f (x0 +h)+f (x0 +h) = 2f (x0 )+h2 f ′′ (x0 )+ (f (ξ1 )+f (4) (ξ2 ))
24
▸ Rearranging,

1
f ′′ (x0 ) ≈ (f (x0 − h) − 2f (x0 ) + f (x0 + h)))
h2
with error,
h2 (4)
∣f (ξ)∣
12
for some ξ in (x0 − h, x0 + h).
▸ Exercise: A five point formula for second derivative.

11 / 33
Round-Off error Instability

▸ Consider the three point central difference formula for first

order derivative:
f (x0 + h) − f (x0 − h) h2 (3)
f ′ (x0 ) = − f (ξ)
2h 6
▸ Now assume that we make following error in roundoff
representations:

f (x0 + h) = f¯(x0 + h) + e(x0 + h)

and
f (x0 − h) = f¯(x0 − h) + e(x0 − h)
where f¯ is the floating point representation of f with error e.

12 / 33
Round-Off error Instability

▸ Then total error in the approximation is,

f (x0 + h) − f (x0 − h) e(x0 + h) − e(x0 − h) h2 (3)

f ′ (x0 )− = − f (ξ)
2h 2h 6
▸ Assuming that the roundoff error e is bounded by ε, we get,

f (x0 + h) − f (x0 − h) ε h2
∣f ′ (x0 ) − ∣ ≤ + M.
2h h 6
▸ What is the optimal value of h?

13 / 33
Round-Off error Instability

▸ Converting it in optimization problem: Find minima of the

function
ε h2
e(h) = + M
h 6
▸ Minima occurs at √
3 3ε
h=
M
▸ Exercise: Repeat the same procedure to get the minimum
value for three point formula for second order derivative.

14 / 33
Topics in this section

1. NewtonCotes Rules
▸ Trapezoidal rule
▸ Simpson’s rule
▸ Composite Rules
References
▸ U. Ascher and C. Greif: Numerical methods. 15.1 − 15.5.
▸ C. Moler: Numerical Computing with Matlab.
▸ W. Gander, M. Gander, and F. Kwok: Scientific
Computing.Springer 2014.

15 / 33
Introduction

Let f (x) be a continuous real-valued function of a real variable,

defined on a finite interval a ≤ x ≤ b. We seek to Compute the
value of the integral,
b
∫ f (x)dx.
a

Quadrature means the approximate evalution of such a definite

integral.
Principle idea:Interpolate f by (piecewise) polynomial, then
integrate the polynomial.

16 / 33
Introduction(Cont.)

The fundamental additive property of definite integrals is the basis

for many quadrature formulae.For a < c < b, we have
b c b
∫ f (x)dx = ∫ f (x)dx + ∫ f (x)dx
a a c

In particular, we can subdivide an interval in many subintervals,

integrate f in the subintervals, and finally add up these partial
integrals to get the overall result.
We can subdivide subintervals that indicate big errors recursively in
an adaptive procedure.

17 / 33
Basic quadrature rules
Let h = b − a be the length of the interval. The two simplest
quadrature rules are
▸ Midpoint rule

a+b
M = hf ( )
2
▸ Trapezoidal rule

f (a) + f (b)
T =h
2
The accuracy of a quadrature rule can be predicted in part by
examining its behavior on polynomials.
Order of a quadrature rule is the degree of the lowest degree
polynomial that the rule does not integrate exactly.
The orders of both midpoint and trapezoidal rule are 2.

18 / 33
Newton-Cotes rule

One obvious way to approximate f is by polynomial interpolation

and then integrate the polynomial.
Let a = x0 < x1 < x2 < ... < xn = b be equidistant points in [a, b].
The Lagrange interpolating polynomial is given by
n n x − xj
Pn (x) = ∑ Li (x)f (xi ), Li (x) = ∏ , i = 0, ..., n.
i=0 j=0 j≠i xi − xj

Then

b b n b n
∫ f (x)dx ≈ ∫ Pn (x)dx = ∑ f (xi ) ∫ Li (x)dx = ∑ f (xi )wi .
a a i=0 a i=0

Clearly, such a method has atleast order n + 1.

19 / 33
Trapezoidal rule again

The trapezoidal rule is the Newton-Cotes rule for n = 1 ∶

b−x x −a
P1 = f (a) + f (b)
b−a b−a
and (with x = a + hξ)

b b−x b x −a 1 h
∫ dx = ∫ dx = h ∫ ξdξ = .
a b−a a b−a 0 2
So,

h
T= (f (a) + f (b)).
2

20 / 33
Simpson’s rule
Simpson’s rule is the Newton-Cotes rule for n = 2 ∶ (x = a + hξ)

P2 (ξ) = f (a)(1 − 2ξ)(1 − ξ) + f (c)4ξ(1 − ξ) + f (b)ξ(2ξ − 1).

h a+b
S = (f (a) + 4f (c) + f (b)), c = .
6 2
Note that

2 1
S= M+ T
3 3
It turns out that Simpson’s rule also integrates cubics exactly, but
not quartics, so its order is four.
Newton-Cotes methods with order up to 8 exist.
Even higher order methods can be computed, but some of their
weights wi get negative which is undesirable for potential round-off.

21 / 33
Summary of basic quadrature rules

Quadrature Rule Error

f ′′ (ξ1 )
Midpoint (b − a)f ( a+b
2 ) 24 (b − a)3

f ′′ (ξ2 )
2 [f (a) + f (b)] − 12 (b − a)3
(b−a)
Trapezoidal

f 4 (ξ3 ) b−a 5
6 [f (a) + 4f ( 2 ) + f (b)] − 90 ( 2 )
(b−a) a+b
Simpson

Exercise: Prove the above error estimates.

22 / 33
Composite Rules

▸ Applying a quadrature rule to ∫ab f (x) may not yield an

approximation with the desired accuracy.
▸ To increase the accuracy, one can partition the interval [a, b]
into subintervals and apply the quadrature rule to each
subinterval.
▸ The resulting formula is known as a composite rule

23 / 33
 Composite Rules(cont.)
mposite Rules (cont.)

Picture from Moler: Numerical Computing with Matlab.

24 / 33
Composite trapezoidal rule

Let [a, b] be partitioned into n equidistant subintervals (xi , xi+1 ) of

length h = xi+1 − xi = (b − a)/n.
Then we apply the trapezoidal rule to each subinterval to obtain
the composite trapezoidal rule
1 1
T (h) = h( y0 + y1 + ... + yn−1 + yn ), yi = f (xi ).
2 2
The error of the composite trapezoidal rule is

RRR b RRR
R (b − a)h ∣ f ′′ (ξ) ∣,
2
RRR
RRR ∫a f (x)dx − T (h)RRRRR = ξ ∈ [a, b].
RR RR 12

Exercise

25 / 33
Composite trapezoidal rule

Let’s assume we have computed T (h). How do we

computeT (h/2)?

h h 1 1
T ( ) = ( y0 + y1 + ... + y2n−1 + y2n ),
2 2 2 2
h 1 1
= ( y0 + y2 + ... + y2n−2 + y2n )
2 2 2
h
+ (y1 + y3 + ... + y2n−1 )
2
1 h
= T (h) + (y1 + y3 + ... + y2n−1 )
2 2

26 / 33
function T=trapez(f,a,b,tol)
Numerical Methods for Computational Science and Engineering
Composite Rules
Trapezoidal rule

function T=trapez(f,a,b,tol);
% TRAPEZ(f,a,b,tol) tries to integrate int_a^b f(x) dx
% to a relative tolerance tol using the composite
% trapezoidal rule.
%
h = b-a; s = (f(a)+f(b))/2;
tnew = h*s; zh = 1; told = 2*tnew;
fprintf(’h = 1/%2.0f, T(h) = %12.10f\n’,1/h, tnew)
while abs(told-tnew)>tol*abs(tnew),
told = tnew; zh = 2*zh; h = h/2;
s = s+sum(f(a+[1:2:zh]*h));
tnew = h*s;
fprintf(’h = 1/%2.0f, T(h) = %12.10f\n’,1/h, tnew)
end;
T = tnew;

NumCSE, Lecture 23, Dec 5, 2013 14/40

27 / 33
Example
1 xe x e −2
I =∫ dx = = 0.3591409142...
0 (x + 1)2 2

Intermediate approximations calculated by trapez(f,0,1,1e-4)

T (hi )−I
i hi T (hi ) T (hi−1 )−I
0 1 0.339785228 .
1
1 2 0.353083866 0.31293
1
2 4 0.357515195 0.26840
1
3 8 0.358726477 0.25492
1
4 16 0.359036783 0.25125
1
5 32 0.359114848 0.25031
1
6 64 0.359134395 0.25007

If f ′′ does not vary too much, error decreases by factor 4/step.

28 / 33
Composite Simpson’s rule

For the Composite Simpson Rule, we need to partition the interval

[a, b] into 2n subintervals.
2n , xi = a + ih and yi = f (xi ) for i = 0, 1, 2, ..., 2n.Applying
Let h = b−a
Simpson’s Rule to two consecutive subintervals
xi+1 h
∫ f (x)dx ≈ (yi−1 + 4yi + yi+1 )
xi−1 3
We obtain the Composite Simpson rule

h
S(h) = (y0 + 4y1 + 2y2 + 4y3 + ... + 2y2n−2 + 4y2n−1 + y2n ).
3

29 / 33
Composite Simpson’s rule(cont.)

The integration error is

RRR b RRR
R (b − a)h ∣f (4) (ξ)∣
4
RRR
RRR ∫a f (x)dx − S(h)RRRRR =
RR RR 180
with a ≤ ξ ≤ b.

Exercise

30 / 33
Composite Simpson’s rule(cont.)

Again Matlab function that computes approximations to the

integral by halving the integration step until a given tolerance is
met.
We introduce variables for sums of function values with the same
weight ∣s1, s2ands4∣ then we have the following update formula to
compute the new Simpson approximation:

S(h) = 3 (s1
h
+ 4s4 + 2s2 )
s1new = s1
s2new = s2 + s4
s4new = sum of new function values
h/2 new
S(h/2) = 3 (s1 + 4s4new + 2s2new )

31 / 33
function S=Simpson(f,a,b,tol);
Numerical Methods for Computational Science and Engineering
Composite Rules
Simpson’s rule

function S=simpson(f,a,b,tol);
% SIMPSON (f,a,b,tol) tries to integrate int_a^b f(x) dx
% to a relative tolerance tol using the composite
% Simpson rule.
%
h = (b-a)/2; s1 = f(a)+f(b); s2 = 0;
s4 = f(a+h); snew = h*(s1+4*s4)/3;
zh = 2; sold = 2*snew;
fprintf(’h = 1/%2.0f, S(h) = %12.10f\n’,1/h, snew)
while abs(sold-snew)>tol*abs(snew),
sold = snew; zh = 2*zh; h = h/2; s2 = s2+s4;
s4 = sum(f(a +[1:2:zh]*h));
snew = h*(s1+2*s2+4*s4)/3;
fprintf(’h = 1/%2.0f, S(h) = %12.10f\n’,1/h, snew)
end
S = snew;

NumCSE, Lecture 23, Dec 5, 2013 19/40

32 / 33
Same example
1 xe x e −2
I =∫ dx = = 0.3591409142...
0 (x + 1) 2 2

Intermediate approximations calculated by Simpson(f , 0, 1, 1e − 8)

S(hi )−I
i hi S(hi ) S(hi−1 )−I
0 1 0.357516745 .
1
1 2 0.358992305 0.09149
1
2 4 0.359130237 0.07184
1
3 8 0.359140219 0.06511
1
4 16 0.359140870 0.06317
1
5 32 0.359140911 0.06267
1
6 64 0.359140914 0.06254

33 / 33