Numerical Methods and Computation

MTL 107

Harish Kumar
(hkumar@iitd.ac.in)

Dept. of Mathematics, IIT Delhi

Introduction

Topics in this Section

▶ Approximation of Derivatives if function values are knows.
▶ Higher order approximation.
▶ Approximation of Second order derivatives.
▶ Richardson Extrapolation
▶ Numerical Differentiation using Lagrange Interpolation
▶ Roundoff Errors and Numerical Differentialtions
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 .
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
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
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 .
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) h3 (3) h3
 
1
f (ξ) = f (ξ2 ) + f (3) (ξ1 )
6 2 3! 3!

by Intermediate Value Theorem.

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)
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.
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.
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.
Richardson Extrapolation

▶ The idea is to generate high order formula from lower order

one.
▶ We will demonstrate this with an example for second order
approximation.
▶ Using Taylor’s theorem and preserving higher order terms we
get,
1
f ′′ (x0 ) = (f (x0 − h) − 2f (x0 ) + f (x0 + h))
h2
h2 (4) h4 (6)
− f (x0 ) − f (x0 ) + O(h5 )
12 360
Richardson Extrapolation

▶ Changing h to 2h we get,

1
f ′′ (x0 ) = (f (x0 − 2h) − 2f (x0 ) + f (x0 + 2h))
4h2
4h2 (4) 16h4 (6)
− f (x0 ) − f (x0 ) + O(h5 ).
12 360
▶ Multiplying first formula with 4, subtracting the other one and
dividing by 3, we get,
1
f ′′ (x0 ) =
12h2
(−f (x0 − 2h) + 16f (x0 − h) − 30f (x0 ) + 16f (x0 + h) − f (x0 + 2h))
Richardson Extrapolation

▶ The error term can be written as,

h4 (4)
|f (ξ)|,
90
for some ξ in (x0 − 2h, x0 + 2h).
▶ Exercise: Repeat the same procedure to derive a five point
central difference formula from three point central difference
formulas for first order derivatives.
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.
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?
Round-Off error Instability

▶ Converting it in optimization problem: Find minima of the

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

1. NewtonCotes Rules
▶ Trapezoidal rule
▶ Simpson’s rule
2. Romberg integration
3. Gauss quadrature
4. Adaptive quadrature
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.
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,
Z 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.
Introduction(Cont.)

The fundamental additive property of definite integrals is the basis

for many quadrature formulae.For a < c < b, we have
Z b Z c Z 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.
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.
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 Y x − xj
Pn (x) = Li (x)f (xi ), Li (x) = , i = 0, ..., n.
xi − xj
i=0 j=0j̸=i

Then

Z b Z b n
X Z b n
X
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.

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 1
b−x x −a
Z Z Z
h
dx = dx = h ξdξ = .
a b−a a b−a 0 2

So,

h
T = (f (a) + f (b)).
2
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.
Summary of basic quadrature rules

Quadrature Rule Error

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

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

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

Exercise: Prove the above error estimates.

Composite Rules

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

R

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
 Composite Rules(cont.)
mposite Rules (cont.)

Picture from Moler: Numerical Computing with Matlab.

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

b
(b − a)h2
Z
f (x)dx − T (h) = | f ′′ (ξ) |, ξ ∈ [a, b].
a 12

Exercise
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
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

Example
1
xe x e −2
Z
I = 2
dx = = 0.3591409142...
0 (x + 1) 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.

Composite Simpson’s rule

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

[a, b] into 2n subintervals.
Let h = b−a2n , xi = a + ih and yi = f (xi ) for
i = 0, 1, 2, ..., 2n.Applying Simpson’s Rule to two consecutive
subintervals
Z 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
Composite Simpson’s rule(cont.)

The integration error is

b
(b − a)h4 (4)
Z
f (x)dx − S(h) = f (ξ)
a 180
with a ≤ ξ ≤ b.

Exercise
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:
h
S(h) = 3 (s1
+ 4s4 + 2s2 )
s1new = s1
s2new = s2 + s4
s4new = sum of new function values
h/2 new
S(h/2) = 3 (s1 + 4s4new + 2s2new )
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

Same example
1
xe x e −2
Z
I = 2
dx = = 0.3591409142...
0 (x + 1) 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
Romberg method

The asymptotic expansion for the approximation of an integral by

the trapezoidal rule is given by
Z b
T (h) = g (t)dt + K1 h2 + K2 h4 + ... + Km h2m + Rm .
a

We can extrapolate to the limit T (0) by interpolating the values

T (hi ) for some hi > 0 by a polynomial.
Error term Rm is small⇒ values T (hi ) are close to the values of
polynomial
Z b
2 4 2m
P2m (h) = K0 + K1 h + K2 h + ... + Km h , with K0 = g (t)dt.
a
Romberg method(cont.)

Only even powers of h ⇒ consider polynomial in variable x = h2 .

We extrapolate to x = 0 the data:
x h02 h12 ... hi2 ...
y T (h0 ) T (h1 ) ... T (hi ) ...

Let Ri,0 = T (hi ) and hi+1 = hi /2. Then, for the errors,

Ei,0 = K0 − Ri,0 = K1 hi2 + K2 hi4 + K3 hi6 + ...

Ei+1,0 = K0 − Ri+1,0 = K1 (hi /2)2 + K2 (hi /2)4 + K3 (hi /2)6 + ...

Eliminate K1 − term :
3 15
Ei,0 − 4Ei+1,0 = K2 hi4 + K3 hi6 + ...
4 16
Romberg method(cont.)

The Corresponding approximations are

4Ri+1,0 − Ri,0
Ri+1,1 =
3
Z b
with Ri+1,1 − g (t)dt = O(h4 )
a

This procedure eliminating higher and higher order error terms can
be continued, yielding the Romberg method:

4j Ri+1,j−1 − Ri,j−1
Ri+1,j =
4j − 1
Romberg method(cont.)

This is the general scheme:

T (h0 ) = R00
T (h1 ) = R10 R11
.. .. .. ..
. . . .
T (hk ) = Rk0 Rk1 ··· Rkk

with
4j R j+1,j−1 − Ri,j−1
Ri+1,j = 0 ≤ j ≤ i ≤ k.
4j − 1
Same example again

1
xe x e −2
Z
I = dx = = 0.3591409142...
0 (x + 1)2 2

Intermediate approximations calculated by

[a z] = romberg (f , 0, 1, 1e − 8)
Same example again(cont.)

Romberg scheme columns 0 to 3

0.33978522855738
0.35308386657870 0.35751674591915
0.35751519587192 0.35899230563633 0.35909067628414
0.35872647716421 0.35913023759497 0.35913943305888 0.35914020697
0.35903678355577 0.35914021901962 0.35914088444793 0.35914090748
0.35911484861929 0.35914087030714 0.35914091372630 0.35914091419
0.35913439576246 0.35914091147685 0.35914091422149 0.35914091422

columns 4 to 6
0.35914091023295
0.35914091421734 0.35914091422123
0.35914091422950 0.35914091422951 0.35914091422952
Another example

R1√ 2
Let’s integrate 0 xdx = 3 by Romberg method

[I T]=romberg(@sqrt,0,1,1e-8)

This generates a warning message that 15 extrapolation steps were

not sufficient to meet the tolerance of 10−8
romberg is hardly more accurate than Simpson!
√
Reason: the function f (x) = x is not differentiable at x = 0 ⇒
the asymptotic expansion T (0) = K0 + K1 h2 + ... does not hold.
Gauss Quadrature

▶ So far: quadrature rules for fixed location of nodes

(ξi = a + ih) at which the function was evaluated.
▶ Gauss quadrature rules allow the nodes to vary as well
▶ The additional freedom can be used to increase the accuracy
of the quadrature rule.
▶ We consider Gauss quadrature rules in the interval [−1, 1]. A
transformation
b−a a+b 2x − b − a
[a, b] ∋ x = t+ ⇐⇒ t = ∈ [−1, 1]
2 2 b−a
may be needed.
b 1  
b−a b−a a+b
Z Z
f (x)dx = f t+ dt.
a 2 −1 2 2
Gauss Quadrature(cont.)

For a given number of nodes n, we want to find the location of

nodes ξk and the weights wk in the quadrature rule
Z 1 n
X
f (x)dx ≈ wk f (ξk ) (Gauss)
−1 k=1

that allows us to integrate exactly a polynomial of degree as high

as possible.
Theorem
The order m of an n-point quadrature rule (Gauss) is bounded by
2n − 1.
Gauss Quadrature(cont)

We consider the case n = 2 for deriving conditions on nodes and

weights.
Four degrees of freedom (2 nodes ξ1 , ξ2 /2 weights w 1, w 2).
We can satisfy 4 nonlinear equations: We derive a quadrature rule
that is exact for polynomial up to degree 3:

Z 1 Z 1
dx = 2 = 1w1 + 1w2 xdx = 0 = ξ1 w1 + ξ2 w2
−1 −1
1 1
2
Z Z
x dx = = ξ12 w1 + ξ22 w2
2
x 3 dx = 0 = ξ13 w1 + ξ23 w2
−1 3 −1
Gauss Quadrature(cont)

The solution is
√ √
{w1 = 1, w2 = 1, ξ1 = 1/ 3, ξ2 = −1/ 3}

and the corresponding Gauss-Legendre rule for n = 2 is

Z 1    
1 1
f (x)dx ≈ f √ + f −√
−1 3 3

This approach(solving a nonlinear system of 2n equations) doesn’t

work for big n.
Gauss quadrature, general case

Remember polynomial interpolation:Determine polynomial pn of

degree at most n such that

f (xi ) = pn (xi ), −1 ≤ x0 < x1 < ...xn ≤ 1.

The basic rule is

Z 1 Z 1
f (x)dx ≈ pn (x)dx.
−1 −1

By means of the interpolation error we have

Z 1 Z 1 n
Y
(f (x) − pn (x))dx = f [x0 , x1 , ..., xn , x] (x − xi )dx
−1 −1 i=1
Gauss quadrature, general case(cont.)
If f ∈ Pm with m ≤ n then f [x0 , x1 , ..., xn , x] = 0 =⇒ rule is exact.
If f ∈ Pm with m > n then f [x0 , x1 , ..., xn , x] ∈ Pm−n−1 .
Best approximation/orthogonal polynomials:Legendre polynomials
Φ1 (x), Φ2 (x), ... are orthogonal w.r.t. inner product
Z 1
f (x)g (x)dx. (0.1)
−1

Then,
Z 1
g (x)Φn+1 (x)dx = 0, ∀g ∈ Pn .
−1

In particular,if f ∈ P2n+1 then

Z 1
f [x0 , x1 , ..., xn , x]Φn+1 (x)dx = 0.
−1
Gauss quadrature,general case(cont.)

Choose zeros x0 ≤ x1 ≤ ... ≤ xn of Legendre polynomial Φn+1

(Gauss points),cf.Chebyshev interpolation.
The weights are given by
1 1 n
x − xj
Z Z Y
wi = Li (x)dx = dx
−1 −1 j=0 xi − xj
j̸=i

where the Li are Lagrange polynomials.

Disadvantage of Gaussian quadrature:integration points xi depend
on n.
▶ Function values f (xi ) cannot be reused if n is increased
▶ Function values f (xi ) cannot be reused if composite rule is
refined.
Adaptive Quadrature

▶ So far: we used composite rules in which all integration steps

had same size h.
▶ Controlled integration error by halving the step size h ← h/2.
▶ We generated equidistant subintervals and applied the same
quadrature rule to each subinterval.
▶ In each subinterval there will be a different error depending on
the local behavior of the function.
▶ The global error will be ‘localized’
▶ To minimize the number of function evaluations, need to
adapted length of subintervals to local behavior of integrand.
▶ The partition should be such that the error becomes roughly
the same for each subinterval.
Adaptive Quadrature:Example
R1√
Let’s again consider 0 xdx = 2/3.
T=trapez(@sqrt,0,1,1e-5) need 2048 function evaluations.
T=trapez(@sqrt,0,1,1e-10) needs 4194304 function evalution
to complete.
Most of the error is in the leftmost interval. For the second
subinterval [h, 2h]
we have,
maxh≤x≤2h f ′′ (x) = h−3/2 /4

f ′′ (ξ)
h3/2
I − T (h) = h3 ≤ = 6.357.10−7 < 10−5 .
12 48

For the subinterval [1 − h, 1], the analogous error estimate gives

the bound 1.94.10−11 .
Adaptive Quadrature:Example(cont)

In an adaptive approach one proceeds as follows (according to the

principle of divide and conquer):
▶ First, integrate f using two different numerical integration
schemes (or the same method with two different step sizes) to
get the approximations I1 and I2 .
▶ f I1 ≈ I2 , one accepts I1 as the value of integral
(|I2 − I1 | < tol)
▶ Otherwise, the interval [a, b] is divided, e.g., into two equal
parts [a, m] and [m, b], where m = (a + b)/2, and the two
respective integrals are computed (conquered) independently.
Adaptive Quadrature:Example(cont)
Numerical Methods for Computational Science and Engineering
Adaptive Quadrature

Adaptive Quadrature: Example (cont.)

function [Q,fcount] = quadtx(F,a,b,tol,varargin)

...

% Initialization
c = (a + b)/2;
fa = F(a,varargin{:});
fc = F(c,varargin{:});
fb = F(b,varargin{:});

% Recursive call
[Q,k] = quadtxstep(F, a, b, tol, fa, fc, fb, varargin{:});
fcount = k + 3;

From Moler: Numerical Computing with Matlab.

NumCSE, Lecture 23, Dec 5, 2013 38/40
Adaptive Quadrature:Example(cont)
Numerical Methods for Computational Science and Engineering
Adaptive Quadrature

Adaptive Quadrature: Example (cont.)

function [Q,fcount] = quadtxstep(F,a,b,tol,fa,fc,fb,varargin)
% Recursive subfunction used by quadtx.
h = b - a;
c = (a + b)/2;
fd = F((a+c)/2,varargin{:});
fe = F((c+b)/2,varargin{:});
Q1 = h/6 * (fa + 4*fc + fb);
Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);
if abs(Q2 - Q1) <= tol
Q = Q2 + (Q2 - Q1)/15;
fcount = 2;
else
[Qa,ka] = quadtxstep(F, a, c, tol, fa, fd, fc, varargin{:});
[Qb,kb] = quadtxstep(F, c, b, tol, fc, fe, fb, varargin{:});
Q = Qa + Qb;
fcount = ka + kb + 2;
end

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

Adaptive Quadrature:Example(cont)

How should tol be chosen in the two subroutines calls in

quadtxstep?
tol/2 in each subinterval? This is too conservative.
Errors in general overestimated. (Extrapolated Simpson’s rule)
We can allow one of the two recursive calls to have an error close
to the tolerance because the other subinterval will probably have a
much smaller error. For these reasons, the same value of tol is
used in each recursive call.
See choice of tol in quadtxstep.