Numerical Analysis and Computing

Composite Quadrature
Romberg Quadrature
Adaptive Quadrature
Gaussian Quadrature
Numerical Analysis and Computing

Lecture Notes #08
Numerical Differentiation and Integration
Composite Numerical Integration; Romberg Integration
Adaptive Quadrature / Gaussian Quadrature
Joe Mahaffy,
[email protected]
Department of Mathematics
Dynamical Systems Group
Computational Sciences Research Center
San Diego State University
San Diego, CA 92182-7720
http://www-rohan.sdsu.edu/jmahaffy
Spring 2010
Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (1/43)
Romberg Quadrature
Adaptive Quadrature
Gaussian Quadrature
Outline
1 Composite Quadrature
Divide and Conquer; Example Simpsons Rule
Generalization
Collecting the Error...
2 Romberg Quadrature
Applying Richardsons Extrapolation
Romberg Quadrature Code Outline
3 Adaptive Quadrature
Introduction
Building the Adaptive CSR Scheme
Example...
Putting it Together...
4 Gaussian Quadrature
Ideas...
2-point Gaussian Quadrature
Higher-Order Gaussian Quadrature Legendre Polynomials
Examples: Gaussian Quadrature in Action; HW#7

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Divide and Conquer with Simpsons Rule I/III
The exact solution:

Z 4
e x dx = e 4 e 0 = 53.59815
0
Simpsons Rule with h = 2

Z 4
2
e x dx (e 0 + 4e 2 + e 4 ) = 56.76958.
0 3
The error is 3.17143 (5.92%).

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Divide and Conquer with Simpsons Rule I/III
The exact solution:

Z 4
e x dx = e 4 e 0 = 53.59815
0
Simpsons Rule with h = 2

Z 4
2
e x dx (e 0 + 4e 2 + e 4 ) = 56.76958.
0 3
The error is 3.17143 (5.92%).
Divide-and-Conquer: Simpsons Rule with h = 1
Z 2 Z 4
1 1
x
e dx+ e x dx (e 0 +4e 1 +e 2 )+ (e 2 +4e 3 +e 4 ) = 53.86385
0 2 3 3
The error is 0.26570. (0.50%) Improvement by a factor of 10!
Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Divide and Conquer with Simpsons Rule II/III
The exact solution:

Z 4
e x dx = e 4 e 0 = 53.59815
0
Divide-and-Conquer: Simpsons Rule with h = 1/2
Z 1 Z 2 Z 3 Z 4
1 1
+ + + e x dx (e 0 + 4e 1/2 + e 1 ) + (e 1 + 4e 3/2 + e 2 )
0 1 2 3 6 6
1 2 1
+ (e + 4e 5/2 + e 3 ) + (e 3 + 4e 7/2 + e 4 ) = 53.61622
6 6
The error has been reduced to 0.01807 (0.034%).

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Divide and Conquer with Simpsons Rule II/III
The exact solution:

Z 4
e x dx = e 4 e 0 = 53.59815
0
Divide-and-Conquer: Simpsons Rule with h = 1/2
Z 1 Z 2 Z 3 Z 4
1 1
+ + + e x dx (e 0 + 4e 1/2 + e 1 ) + (e 1 + 4e 3/2 + e 2 )
0 1 2 3 6 6
1 2 1
+ (e + 4e 5/2 + e 3 ) + (e 3 + 4e 7/2 + e 4 ) = 53.61622
6 6
The error has been reduced to 0.01807 (0.034%).
h abs-error err/h err/h2 err/h3 err/h4
2 3.17143 1.585715 0.792857 0.396429 0.198214
1 0.26570 0.265700 0.265700 0.265700 0.265700
1/2 0.01807 0.036140 0.072280 0.144560 0.289120
Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Divide and Conquer with Simpsons Rule III/III
Extending the table...

h abs-error err/h err/h2 err/h3 err/h4 err/h5
2 3.171433 1.585716 0.792858 0.396429 0.198215 0.099107
1 0.265696 0.265696 0.265696 0.265696 0.265696 0.265696
1/2 0.018071 0.036142 0.072283 0.144566 0.289132 0.578264
1/4 0.001155 0.004618 0.018473 0.073892 0.295566 1.182266
1/8 0.000073 0.000580 0.004644 0.037152 0.297215 2.377716
Clearly, the err/h4 column seems to converge (to a non-zero

constant) as h 0. The columns to the left seem to converge to
zero, and the err/h5 column seems to grow.
This is numerical evidence
that the composite Simpsons rule has
a convergence rate of O h4 . But, isnt Simpsons rule 5th
order???

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Generalized Composite Simpsons Rule I/II
For an even integer n: Subdivide the interval [a, b] into n

subintervals, and apply Simpsons rule on each consecutive pair of
sub-intervals. With h = (b a)/n and xj = a + jh, j = 0, 1, . . . , n,
we have
Z b n/2 Z x2j
X
f (x)dx = f (x)dx
a j=1 x2j2
n/2
X h h5
= f (x2j2 ) + 4f (x2j1 ) + f (x2j ) f (4) (j ) ,
3 90
j=1
for some j [x2j2 , x2j ], if f C 4 [a, b].

Since all the interior even x2j points appear twice in the sum, we
can simplify the expression a bit...
Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Generalized Composite Simpsons Rule II/II

Z b n/2
h X
f (x)dx = f (x0 ) f (xn ) + 4f (x2j1 ) + 2f (x2j )
a 3
j=1
n/2
h5 X
f (4) (j ).
90
j=1
The error term is:

n/2
h5 X (4)
E (f ) = f (j ), j [x2j2 , x2j ]
90
j=1

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
The Error for Composite Simpsons Rule I/II
If f C 4 [a, b], the Extreme Value Theorem implies that f (4)

assumes its max and min in [a, b]. Now, since
min f (4) (x) f (4) (j ) max f (4) (x),
x[a,b] x[a,b]

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature

min f (4) (x) f (4) (j ) max f (4) (x),
x[a,b] x[a,b]
hni n/2 hni

X
(4)
min f (x) f (4) (j ) max f (4) (x),
2 x[a,b] 2 x[a,b]
j=1

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature

min f (4) (x) f (4) (j ) max f (4) (x),
x[a,b] x[a,b]
hni n/2 hni

X
(4)
min f (x) f (4) (j ) max f (4) (x),
2 x[a,b] 2 x[a,b]
j=1
n/2
X
(4) 2
min f (x) f (4) (j ) max f (4) (x),
x[a,b] n x[a,b]
j=1

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature

min f (4) (x) f (4) (j ) max f (4) (x),
x[a,b] x[a,b]
hni n/2 hni

X
(4)
min f (x) f (4) (j ) max f (4) (x),
2 x[a,b] 2 x[a,b]
j=1
n/2
X
(4) 2
min f (x) f (4) (j ) max f (4) (x),
x[a,b] n x[a,b]
j=1
By the Intermediate Value Theorem (a, b) so that
n/2 n/2
2 X (4) n (4) X
f (4) () = f (j ) f () = f (4) (j )
n 2
j=1 j=1
Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
The Error for Composite Simpsons Rule II/II
We can now rewrite the error term:

n/2
h5 X (4) h5 (4)
E (f ) = f (j ) = nf (),
90 180
j=1
or, since h = (b a)/n n = (b a)/h, we can write
(b a) 4 (4)
E (f ) = h f ().
180
Hence Composite Simpsons Rule has degree of accuracy 3
(since it is exact for polynomials up to order 3),
and the error is
proportional to h4 Convergence Rate O h4 .

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Composite Simpsons Rule Algorithm
Algorithm (Composite Simpsons Rule)

Given the end points a and b and an even positive integer n:
[1] h = (b a)/n
[2] ENDPTS = f(a)+f(b)
ODDPTS = 0
EVENPTS = 0
[3] FOR i = 1, . . . , n 1 (interior points)
x =a+i h
if i is even: EVENPTS += f(x)
if i is odd: ODDPTS += f(x)
END
[4] INTAPPROX = h*(ENDPTS+2*EVENPTS+4*ODDPTS)/3

Romberg Quadrature
Generalization
Adaptive Quadrature
Gaussian Quadrature
Homework #7 Due Friday 11/13/2009, 12-noon
(Part-1)
Implement Composite Simpsons Rule, and use your code to solve
BF-4.4.3-a,b,c,d.

Romberg Quadrature Applying Richardsons Extrapolation
Adaptive Quadrature Romberg Quadrature Code Outline
Gaussian Quadrature
Romberg Integration The Return of Richardsons Extrapolation
Romberg Integration is the combination of the Composite

Trapezoidal Rule (CTR)

Z b n1
h X (b a) 2
f (x)dx = f (a) + f (b) + 2 f (xj ) h f ()
a 2 12
j=1
and Richardson Extrapolation.

Here, we know that the error term for regular Trapezoidal Rule is
O(h3 ). By the same argument as for Composite Simpsons Rule,
this gets reduced to O(h2 ) for the composite version.

Gaussian Quadrature
Romberg Integration Step-1: CTR Refinement
Let Rk,1 denote the Composite Trapezoidal Rule with 2k1

sub-intervals, and hk = (b a)/2k1 . We get:
h1
R1,1 = 2 [f (a) + f (b)]
h2
R2,1 = 2 [f (a) + 2f (a + h2 ) + f (b)]
(ba)
= 4 [f (a) + f (b) + 2f (a + h2 )]
1
= 2 [R1,1 + h1 f (a + h2 )]
..
.
k2
2X
1
Rk,1 = Rk1,1 + hk1 f (a + (2i 1)hk )
2
i=1
| {z }
Update formula, using previous value + new points

Gaussian Quadrature
R
Example: Rk,1 for 0 sin(x)dx
1 1 1
f(x) f(x) f(x)

R[2,1] R[3,1] R[4,1]
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
0 1 2 3 0 1 2 3 0 1 2 3
k Rk,1
1 0
2 1.5707963267949
3 1.8961188979370
4 1.9742316019455
5 1.9935703437723
6 1.9983933609701
7 1.9995983886400

Gaussian Quadrature
Extrapolate using Richardson
We know that the error term is O(h2 ), so in order to eliminate this

term we combine to consecutive entries Rk1,1 and Rk,1 to form a
higher order approximation Rk,2 of the integral.
Rk,1 Rk1,1
Rk,2 = Rk,1 +
22 1
Rk,1 O h2
`
Rk,2
0 0
1.5707963267949 2.09439510239
1.8961188979370 2.00455975498
1.9742316019455 2.00026916994
1.9935703437723 2.00001659104
1.9983933609701 2.00000103336
1.9995983886400 2.00000006453

Gaussian Quadrature
Extrapolate, again...
It turns out (Taylor expand to check) that the complete error term
for the Trapezoidal rule only has even powers of h:
Z b
X
f (x) = Rk,1 E2i hk2i .
a i=1
Hence the Rk,2 approximations have error terms that are of size
O(h4 ).
To get O(h6 ) approximations, we compute
Rk,2 Rk1,2
Rk,3 = Rk,2 +
42 1

Gaussian Quadrature
Extrapolate, yet again...
In general, since we only have even powers of h in the error

expansion:
Rk,j1 Rk1,j1
Rk,j = Rk,j1 +
4j1 1
R
Revisiting 0 sin(x)dx:
Rk,1 O h2 Rk,2 O h4 Rk,3 O h6 Rk,4 O h8
` ` ` `
0
1.570796326794897 2.094395102393195
1.896118897937040 2.004559754984421 1.998570731823836
1.974231601945551 2.000269169948388 1.999983130945986 2.000005549979671
1.993570343772340 2.000016591047935 1.999999752454572 2.000000016288042
1.998393360970145 2.000001033369413 1.999999996190845 2.000000000059674
1.999598388640037 2.000000064530001 1.999999999940707 2.000000000000229

Gaussian Quadrature
Homework? No, enough already Heres the code outline!
Code (Romberg Quadrature)

% Romberg Integration for sin(x) over [0,pi]
a = 0; b = pi; % The Endpoints
R = zeros(7,7);
R(1,1) = (b a)/2 (sin(a) + sin(b));
for k = 2 : 7
h = (b a)/2(k1) ;
R(k,1)=1/2 (R(k 1, 1) + 2 h (sin(a + (2 (1 : (2(k2) )) 1) h)));
P
end
for j = 2 : 7
for k = j : 7
R(k,j) = R(k, j 1) + (R(k, j 1) R(k 1, j 1))/(4(j1) 1);
end
end
disp(R)

Composite Quadrature Introduction
Romberg Quadrature Building the Adaptive CSR Scheme
Adaptive Quadrature Example...
Gaussian Quadrature Putting it Together...
More Advanced Numerical Integration Ideas
Adaptive and Gaussian Quadrature

Introduction Adaptive Quadrature
The composite formulas require equally spaced nodes.

This is not good if the function we are trying to integrate has both
regions with large fluctuations, and regions with small variations.
1.5
0.5
-0.5
-1
-1.5
-10 0 10
We need many points where the function fluctuates, but few

points where it is close to constant or linear.
Introduction Adaptive Quadrature Methods
Idea Cleverly predict (or measure) the amount of variation and au-
tomatically add more points where needed.
We are going to discuss this in the context of Composite
Simpsons rule, but the approach can be adopted for other
integration schemes.
First we are going to develop a way to measure the error
a numerical estimate of the actual error in the numerical
integration. Note: just knowing the structure of the error
term is not enough! (We will however use the structure of the
error term in our derivation of the numerical error estimate.)
Then we will use the error estimate to decide whether to accept the
value from CSR, or if we need to refine further (recompute
with smaller h).

Some Notation One-step Simpsons Rule S(f ; a, b)
Notation One-step Simpsons Rule:

Z b
h15 (4)
f (x) dx = S(f ; a, b) f (1 ), 1 (a, b),
a |90 {z }
E(f; h1 ,1 )
where

(b a) a+b (b a)
S(f ; a, b) = f (a) + 4f 2 + f (b) , h1 = .
6 2

Composite Simpsons Rule (CSR)
With this notation, we can write CSR with n = 4, and

h2 = (b a)/4 = h1 /2:
Z b
f (x) dx = S(f ; a, a+b a+b
2 ) + S(f ; 2 , b) E (f ; h2 , 2 ).
a
We can squeeze out an estimate for the error by noticing that

1 h15 (4) 1
E (f ; h2 , 2 ) = f (2 ) = E (f ; h1 , 2 ).
16 90 16
Now, assuming f (4) (1 ) f (4) (2 ), we do a little bit of algebra

magic with our two approximations to the integral...

Wait! Wait! Wait! I pulled a fast one!

1 h15 (4) 1 1 h15 (4) 2
E (f ; h2 , 2 ) = f (2 ) + f (2 )
32 90 32 90
where 12 [a, a+b 2 a+b
2 ], 2 [ 2 , b].
If f C 4 [a, b], then we can use our old friend, the intermediate
value theorem:
f (4) (12 ) + f (4) (22 )
2 [12 , 22 ] [a, b] : f (4) (2 ) = .
2
So it follows that

1 h15 (4)
E (f ; h2 , 2 ) = f (2 ) .
16 90

Back to the Error Estimate...
Now we have

1 h15 (4)
S(f ; a, a+b
2 ) + S(f ;
a+b
2 , b) f (2 )
16 90
h15 (4)
= S(f ; a, b) f (1 ).
90
Now use the assumption f (4) (1 ) f (4) (2 ) (and replace 1 and
2 by ):

h51 (4) 16
f () S(f ; a, b)S(f ; a, (a+b)/2)S(f ; (a+b)/2, b) ,
90 15
h15 (4)
notice that 90 f () = E (f ; h1 , ) = 16E (f ; h2 , ). Hence

1
E (f ; h2 , ) S(f ; a, b)S(f ; a, (a+b)/2)S(f ; (a+b)/2, b) ,
15
Finally, we have the error estimate in hand...
h15 (4)
Using the estimate of 90 f (), we have
Error Estimate for CSR

Z b

f(x)dx S(f; a, (a + b)/2) S(f; (a + b)/2, b)
a
1
S(f; a, b) S(f; a, (a + b)/2) S(f; (a + b)/2, b)
15
Notice!!! S(f ; a, (a + b)/2) + S(f ; (a + b)/2, b) approximates

Rb
a f (x)dx 15 times better than it agrees with the
known quantity S(f ; a, b)!!!

Example Error Estimates I/II
We will apply Simpsons rule to

Z /2
sin(x) dx = 1.
0
Here,
S1 (sin(x); 0, /2) = S(sin(x); 0, /2)

= sin(0) + 4 sin(/4) + sin(/2) = 2 2+1
12 12
= 1.00227987749221.
S2 (sin(x); 0, /2) = S(sin(x); 0, /4) + S(sin(x); /4, /2)

= sin(0) + 4 sin(/8) + 2 sin(/4) + 4 sin(3/8) + sin(/2)
24
= 1.00013458497419.
Example Error Estimates II/II
The error estimate is given by

1
S1 (sin(x); 0, /2) S2 (sin(x); 0, /2)
15
1
= 1.00227987749221 1.00013458497419
15
= 0.00014301950120.
This is a very good approximation of the actual error, which is

0.00013458497419.
OK, we know how to get an error estimate. How do we use
this to create an adaptive integration scheme???

Adaptive Quadrature
Rb
We want to approximate I = a f (x) dx with an error less than
(a specified tolerance).
[1] Compute the two approximations
S1 (f (x); a, b) = S(f (x); a, b), and
S2 (f (x); a, b) = S(f (x); a, a+b
2 ) + S(f (x);
a+b
2 , b).
[2] Estimate the error, if the estimate is less than , we are done.
Otherwise...
[3] Apply steps [1] and [2] recursively to the intervals
[a, a+b a+b
2 ] and [ 2 , b] with tolerance /2.

Adaptive Quadrature, Interval Refinement Example #1
The funny figure above is supposed to illustrate a possible

sub-interval refinement hierarchy. Red dashed lines illustrate
failure to satisfy the tolerance, and black lines illustrate satisfied
tolerance.
level tol interval
1 [a, b]
2 /2 [a, a + ba
2
] [a + (b a)/2, b]
ba
3 /4 [a, a + 4
] [a + ba
4
,a + ba
2
]
.
.
.

Adaptive Quadrature, Interval Refinement Example #2
Adaptive CSR The Function Adaptive CSR Refinement Levels

1 20
0.9 18
0.8 16
0.7 14
Refinement Level
0.6 12
f(x)
0.5 10
0.4 8
0.3 6
0.2 4
0.1 2
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x x
q
3 2
Figure: Application of adaptive CSR to the function f (x) = 1 (x 2e
) . Here,
we have required that the estimated error be less than 106 . The left panel shows
the function, and the right panel shows the number of refinement levels needed to
reach the desired accuracy. At completion we have the value of the integral being
0.61692712, with an estimated error of 3.93 107 .

Composite Quadrature Ideas...
Romberg Quadrature 2-point Gaussian Quadrature
Adaptive Quadrature Higher-Order Gaussian Quadrature Legendre Polynomials
Gaussian Quadrature Examples: Gaussian Quadrature in Action; HW#7
Gaussian Quadrature
Idea: Evaluate the function at a set of optimally chosen points

in the interval.
We will choose {x0 , x1 , . . . , xn } [a, b] and coefficients ci , so that

the approximation
Z b n
X
f (x)dx ci f (xi )
a i=0
is exact for the largest class of polynomials possible.

We have already seen that the open Newton-Cotes formulas
sometimes give us better bang-for-buck than the closed formulas
(e.g. the mid-point formula uses only 1 point and is as accurate as
the two-point trapezoidal rule). Gaussian quadrature takes this
one step further.
Quadrature Types A Comparison
Newton-Cotes Gaussian
Open Closed
Quadrature Degree of Degree of Degree of
Points Accuracy Accuracy Accuracy
1 1 1
2 1 1 3
3 3 3# 5
4 3 3 7
5 5 5 9
The mid-point rule.
Trapezoidal rule.
# Simpsons rule.
The mid-point rule is the only optimal scheme we have see so far.

Gaussian Quadrature Example 2-Point Formula
Suppose we want to find an optimal two-point formula:

Z 1
f (x) dx = c1 f (x1 ) + c2 f (x2 ).
1
Since we have 4 parameters to play with, we can generate a
formula that is exact up to polynomials of degree 3. We get
the following 4 equations:
Z 1

1 dx = 2 = c1 + c2
c1 = 1

Z1
1
c2 = 1
x dx = 0 = c1 x1 + c2 x2

Z1
1
3
x1 =
x 2 dx = 2
3 = c1 x12 + c2 x22

3

Z1
1 3
x2 =
x 3 dx = 0 = c1 x13 + c2 x23

3
1
Higher Order Gaussian Quadrature Formulas
We could obtain higher order formulas by adding more points,

computing the integrals, and solving the resulting non-linear
system of equations... but it gets very painful, very fast.
The Legendre Polynomials come to our rescue!
The Legendre polynomials Pn (x) are orthogonal on [1, 1] with
respect to the weight function w (x) = 1, i.e.
Z 1
0 m=6 n
Pn (x)Pm (x) dx = n n,m =
1 n m = n.
If P(x) is a polynomial of degree less than n, then

Z 1
Pn (x)P(x) dx = 0.
1

A Quick Note on Legendre Polynomials
We will see Legendre polynomials in more detail later. For now,

all we need to know is that they satisfy the property
Z 1
Pn (x)Pm (x) dx = n n,m .
1
and the first few Legendre polynomials are
P0 (x) = 1
P1 (x) = x
P2 (x) = x 2 1/3
P3 (x) = x 3 3x/5
P4 (x) = x 4 6x 2 /7 + 3/35
P5 (x) = x 5 10x 3 /9 + 5x/21.
It turns out that the roots of the Legendre polynomials are the
nodes in Gaussian quadrature.
Higher Order Gaussian Quadrature Formulas
Theorem
Suppose that {x1 , x2 , . . . , xn } are the roots of the nth Legendre
polynomial Pn (x) and that for each i = 1, 2, . . . , n, the coefficients
ci are defined by
Z 1 Y n
x xj
ci = dx.
1 xi x j
j =1
j 6= i
If P(x) is any polynomial of degree less than 2n, then

Z 1 n
X
P(x) dx = ci P(xi ).
1 i=1

Proof of the Theorem I/III
Let us first consider a polynomial, P(x) with degree less than n.

P(x) can be rewritten as an (n 1)-st Lagrange polynomial with
nodes at the roots of the nth Legendre polynomial Pn (x). This
representation is exact since the error term involves the nth
derivative of P(x), which is zero. Hence,
Z 1 n Y
X n
Z 1
x xj
P(x) dx = P(xj ) dx
1 1 xi x j
i=1 j = 1
j 6= i
n Z 1 n n
X Y x xj X
= dx P(xj ) = ci P(xi ),
1 xi xj
i=1 j =1 i=1
j 6= i
which verifies the result for polynomials of degree less than n.

Proof of the Theorem II/III
If the polynomial P(x) of degree [n, 2n) is divided by the nth

Legendre polynomial Pn (x), we get:
P(x) = Q(x)Pn (x) + R(x)
where both Q(x) and R(x) are of degree less than n.

[1] Since deg(Q(x)) < n
Z 1
Q(x)Pn (x) dx = 0.
1
[2] Further, since xi is a root of Pn (x):
P(xi ) = Q(xi )Pn (xi ) + R(xi ) = R(xi ).

Proof of the Theorem III/III
[3] Now, since deg(R(x)) < n, the first part of the proof implies
Z 1 X n
R(x) dx = ci R(xi ).
1 i=1
Putting [1], [2] and [3] together we arrive at
Z 1 Z 1

P(x) dx = Q(x)Pn (x) + R(x) dx
1 1
Z 1 Xn
= R(x) dx = ci R(xi )
1 i=1
n
X
= ci P(xi ),
i=1
which shows that the formula is exact for all polynomials P(x) of
degree less than 2n.
Gaussian Quadrature beyond the interval [1, 1]
By a simple linear transformation,
2x a b (b a)t + (b + a)
t= x= ,
ba 2
we can apply the Gaussian Quadrature formulas to any interval
Z b Z 1
(b a)t + (b + a) (b a)
f (x) dx = f dt.
a 1 2 2

Examples I/II
Degree Pn (x) Roots / Quadrature

points

2 x 2 1/3 {1/ 3, p1/ 3}
3
p
3 x 3x/5 { 3/5, 0, 3/5}
4 x 4 6x 2 /7 + 3/35 {0.86114, 0.33998, 0.33998, 0.86114}
Table: Quadrature points on standard interval.
Z /4
1
(cos(x))2 dx = + = 0.642699081698724
0 4 8
Degree Quadrature points Coefficients

2 0.16597, 0.61942 1, 1
3 0.08851, 0.39270, 0.69688 0.55556, 0.88889, 0.55556
4 0.05453, 0.25919, 0.52621, 0.73087 0.34785, 0.65215, 0.65215, 0.34785
Table: Quadrature points translated to interval of interest; with weight coefficients.

Examples II/II
Z /4
1
(cos(x))2 dx = + = 0.642699081698724
0 4 8
Degree Quadrature points Coefficients

2 0.16597, 0.61942 1, 1
3 0.08851, 0.39270, 0.69688 0.55556, 0.88889, 0.55556
4 0.05453, 0.25919, 0.52621, 0.73087 0.34785, 0.65215, 0.65215, 0.34785
Table: Quadrature points translated to interval of interest; with weight coefficients.
Degree Integral approximation Error

2 0.642317235049753 0.0003818466489...
3 0.642701112090729 0.0000020303920...
4 0.642699075999924 0.0000000056988...
Table: Approximation and Error, for GQ.

Numerical Analysis and Computing

Uploaded by

Copyright:

Available Formats

Numerical Analysis and Computing

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Analysis and Computing

Uploaded by

Copyright:

Available Formats

Composite Quadrature

Numerical Analysis and Computing

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (2/43)

Divide and Conquer with Simpsons Rule I/III

The exact solution:

Simpsons Rule with h = 2

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (3/43)

Divide and Conquer with Simpsons Rule I/III

The exact solution:

Simpsons Rule with h = 2

Divide and Conquer with Simpsons Rule II/III

The exact solution:

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (4/43)

Divide and Conquer with Simpsons Rule II/III

The exact solution:

Divide and Conquer with Simpsons Rule III/III

Extending the table...

Clearly, the err/h4 column seems to converge (to a non-zero

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (5/43)

Generalized Composite Simpsons Rule I/II

For an even integer n: Subdivide the interval [a, b] into n

for some j [x2j2 , x2j ], if f C 4 [a, b].

Generalized Composite Simpsons Rule II/II

The error term is:

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (7/43)

The Error for Composite Simpsons Rule I/II

If f C 4 [a, b], the Extreme Value Theorem implies that f (4)

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (8/43)

The Error for Composite Simpsons Rule I/II

If f C 4 [a, b], the Extreme Value Theorem implies that f (4)

hni n/2 hni

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (8/43)

The Error for Composite Simpsons Rule I/II

If f C 4 [a, b], the Extreme Value Theorem implies that f (4)

hni n/2 hni

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (8/43)

The Error for Composite Simpsons Rule I/II

If f C 4 [a, b], the Extreme Value Theorem implies that f (4)

hni n/2 hni

The Error for Composite Simpsons Rule II/II

We can now rewrite the error term:

or, since h = (b a)/n n = (b a)/h, we can write

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (9/43)

Composite Simpsons Rule Algorithm

Algorithm (Composite Simpsons Rule)

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (10/43)

Homework #7 Due Friday 11/13/2009, 12-noon

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (11/43)

Romberg Integration The Return of Richardsons Extrapolation

Romberg Integration is the combination of the Composite

and Richardson Extrapolation.

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (12/43)

Romberg Integration Step-1: CTR Refinement

Let Rk,1 denote the Composite Trapezoidal Rule with 2k1

Joe Mahaffy, [email protected] Composite-Romberg-Adaptive-Gaussian (13/43)

f(x) f(x) f(x)

0.6 0.6 0.6

0.4 0.4 0.4