Chapter 4
Chapter 4
Chapter 4
𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 )
𝑓 ′ (𝑥0 ) = lim ,
ℎ→0 ℎ
choose a small ℎ and approximate
𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 )
𝑓 ′ (𝑥0 ) ≈
ℎ
The error term for the linear Lagrange polynomial gives:
𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 ) ℎ ″
𝑓 ′ (𝑥0 ) = − 𝑓 (𝜉)
ℎ 2
Also known as the forward-difference formula if ℎ > 0 and the
backward-difference formula if ℎ < 0
General Derivative Approximations
to get
𝑛
′
𝑓 (𝑛+1) (𝜉(𝑥𝑗 ))
𝑓 (𝑥𝑗 ) = ∑ 𝑓(𝑥𝑘 )𝐿′𝑘 (𝑥𝑗 ) + ∏(𝑥 − 𝑥𝑘 )
𝑘=0
(𝑛 + 1)! 𝑘≠𝑗 𝑗
1 ℎ2
𝑓 ′ (𝑥0 ) = [𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 − ℎ)] − 𝑓 (3) (𝜉1 )
2ℎ 6
Suppose that round-off errors 𝜀 are introduced when computing
𝑓. Then the approximation error is
̃ + ℎ) − 𝑓(𝑥
𝑓(𝑥 ̃ − ℎ) 𝜀 ℎ2
0 0
∣𝑓 ′ (𝑥0 ) − ∣ ≤ + 𝑀 = 𝑒(ℎ)
2ℎ ℎ 6
𝑀 − 𝑁 (ℎ) = 𝐾1 ℎ + 𝐾2 ℎ2 + 𝐾3 ℎ3 + ⋯
𝑀 − 𝑁 (ℎ) = 𝐾2 ℎ2 + 𝐾4 ℎ4 + ⋯
Simpson’s Rule
Second Lagrange polynomial with 𝑥0 = 𝑎, 𝑥2 = 𝑏, 𝑥1 = 𝑎 + ℎ,
ℎ = (𝑏 − 𝑎)/2 gives
𝑥2
ℎ ℎ5
∫ 𝑑𝑥 = [𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )] − 𝑓 (4) (𝜉)
𝑥0
3 90
Definition
The degree of accuracy, or precision, of a quadrature formula is the
largest positive integer 𝑛 such that the formula is exact for 𝑥𝑘 , for
each 𝑘 = 0, 1, … , 𝑛.
The Newton-Cotes Formulas
The Closed Newton-Cotes Formulas
Use nodes 𝑥𝑖 = 𝑥0 + 𝑖ℎ, 𝑥0 = 𝑎, 𝑥𝑛 = 𝑏, ℎ = (𝑏 − 𝑎)/𝑛:
𝑏 𝑛
∫ 𝑓(𝑥) 𝑑𝑥 ≈ ∑ 𝑎𝑖 𝑓(𝑥𝑖 )
𝑎 𝑖=0
𝑥𝑛 𝑥𝑛 (𝑥 − 𝑥𝑗 )
𝑎𝑖 = ∫ 𝐿𝑖 (𝑥) 𝑑𝑥 = ∫ ∏ 𝑑𝑥
𝑥0 𝑥0 𝑗≠𝑖
(𝑥𝑖 − 𝑥𝑗 )
Theorem
Let 𝑓 ∈ 𝐶 4 [𝑎, 𝑏], 𝑛 even, ℎ = (𝑏 − 𝑎)/𝑛, 𝑥𝑗 = 𝑎 + 𝑗ℎ, 𝜇 ∈ (𝑎, 𝑏).
The Composite Simpson’s rule for 𝑛 subintervals is
𝑏 (𝑛/2)−1 𝑛/2
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑓(𝑎) + 2 ∑ 𝑓(𝑥2𝑗 ) + 4 ∑ 𝑓(𝑥2𝑗−1 ) + 𝑓(𝑏)]
𝑎
3 𝑗=1 𝑗=1
𝑏 − 𝑎 4 (4)
− ℎ 𝑓 (𝜇)
180
Romberg Integration
ℎ1 (𝑏 − 𝑎)
𝑅1,1 = [𝑓(𝑎) + 𝑓(𝑏)] = [𝑓(𝑎) + 𝑓(𝑏)]
2 2
1
𝑅2,1 = [𝑅1,1 + ℎ1 𝑓(𝑎 + ℎ2 )]
2
1
𝑅3,1 = {𝑅2,1 + ℎ2 [𝑓(𝑎 + ℎ3 ) + 𝑓(𝑎 + 3ℎ3 )]}
2
2𝑘−2
1
𝑅𝑘,1 = [𝑅𝑘−1,1 + ℎ𝑘−1 ∑ 𝑓(𝑎 + (2𝑖 − 1)ℎ𝑘 )]
2 𝑖=1
𝑅𝑘,𝑗−1 − 𝑅𝑘−1,𝑗−1
𝑅𝑘,𝑗 = 𝑅𝑘,𝑗−1 +
4𝑗−1 − 1
Romberg Integration
MATLAB Implementation
function R = romberg(f, a, b, n)
% Compute integral of f(x) from a to b using Romberg integration.
h = b-a;
R = zeros(n,n);
R(1,1) = h/2 * (f(a) + f(b));
for i = 2:n
R(i,1) = 1/2 * (R(i-1,1) + h*sum(f(a + ((1:2^(i-2))-0.5)*h)));
for j = 2:i
R(i,j) = R(i,j-1) + (R(i,j-1)-R(i-1,j-1)) / (4^(j-1)-1);
end
h = h/2;
end
Error Estimation
The assumption 𝑓 (4) (𝜇) ≈ 𝑓 (4) (𝜇)̃ gives the error estimate
𝑏
𝑎+𝑏 𝑎+𝑏
∣∫ 𝑓(𝑥) 𝑑𝑥 − 𝑆 (𝑎, ) − 𝑆( , 𝑏)∣
𝑎
2 2
1 𝑎+𝑏 𝑎+𝑏
≈ ∣𝑆(𝑎, 𝑏) − 𝑆 (𝑎, ) − 𝑆( , 𝑏)∣
15 2 2
Adaptive Quadrature
𝑃0 (𝑥) = 1, 𝑃1 (𝑥) = 𝑥
1 3
𝑃2 (𝑥) = 𝑥2 − 𝑃3 (𝑥) = 𝑥3 − 𝑥
3 5
6 3
𝑃4 (𝑥) = 𝑥4 − 𝑥2 +
7 35
Gaussian Quadrature
Theorem
Suppose 𝑥1 , … , 𝑥𝑛 are roots of 𝑃𝑛 (𝑥) and
1 𝑛 𝑥 − 𝑥𝑗
𝑐𝑖 = ∫ ∏ 𝑑𝑥
−1 𝑗≠𝑖
𝑥𝑖 − 𝑥𝑗
MATLAB Implementation
function [x, c] = gaussquad(n)
% Compute Gaussian quadrature points and coefficients.
P = zeros(n+1,n+1);
P([1,2],1) = 1;
for k = 1:n-1
P(k+2,1:k+2) = ((2*k+1)*[P(k+1,1:k+1) 0] - ...
k*[0 0 P(k,1:k)]) / (k+1);
end
x = sort(roots(P(n+1,1:n+1)));
A = zeros(n,n);
for i = 1:n
A(i,:) = polyval(P(i,1:i),x)';
end
c = A \ [2; zeros(n-1,1)];
Arbitrary Intervals
2𝑥 − 𝑎 − 𝑏 1
𝑡= ⇔ 𝑥 = [(𝑏 − 𝑎)𝑡 + 𝑎 + 𝑏]
𝑏−𝑎 2
Gaussian quadrature then gives
𝑏 1
(𝑏 − 𝑎)𝑡 + (𝑏 + 𝑎) (𝑏 − 𝑎)
∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑓 ( ) 𝑑𝑡
𝑎 −1
2 2
Double Integrals
(𝑑 − 𝑐)(𝑏 − 𝑎) 4 𝜕 4 𝑓 4
4𝜕 𝑓
𝐸=− [ℎ ( 𝜂,
̄ 𝜇)
̄ + 𝑘 (𝜂,̂ 𝜇)]
̂
180 𝜕𝑥4 𝜕𝑦4
𝑑 •1 •4 •2 •4 •1
•4 •16 •8 •16 •4
𝑐 •1 •4 •2 •4 •1
𝑎 𝑏
Non-Rectangular Regions
The step size for 𝑥 is still ℎ = (𝑏 − 𝑎)/𝑛, but for 𝑦 it varies with 𝑥:
𝑑(𝑥) − 𝑐(𝑥)
𝑘(𝑥) =
𝑚
Gaussian Double Integration
More generally, if
𝑔(𝑥)
𝑓(𝑥) = , 0 < 𝑝 < 1, 𝑔 continuous on [𝑎, 𝑏],
(𝑥 − 𝑎)𝑝
𝑔″ (𝑎)
𝑃4 (𝑥) = 𝑔(𝑎) + 𝑔′ (𝑎)(𝑥 − 𝑎) + (𝑥 − 𝑎)2
2!
𝑔‴ (𝑎) 𝑔(4) (𝑎)
+ (𝑥 − 𝑎)3 + (𝑥 − 𝑎)4
3! 4!
Improper Integrals with a Singularity
and write
𝑏 𝑏 𝑏
𝑔(𝑥) − 𝑃4 (𝑥) 𝑃4 (𝑥)
∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑝
𝑑𝑥 + ∫ 𝑑𝑥
𝑎 𝑎
(𝑥 − 𝑎) 𝑎
(𝑥 − 𝑎)𝑝
For the first integral, use the Composite Simpson’s rule to compute
the integral of 𝐺 on [𝑎, 𝑏], where
𝑔(𝑥)−𝑃4 (𝑥)
(𝑥−𝑎)𝑝 , if 𝑎 < 𝑥 ≤ 𝑏
𝐺(𝑥) = {
0, if 𝑥 = 𝑎
Note that 0 < 𝑝 < 1 and 𝑃4 (𝑎) agrees with 𝑔(𝑘) (𝑎) for each
(𝑘)
which gives
∞ 0 1/𝑎
1 𝑡𝑝 1
∫ 𝑝
𝑑𝑥 = ∫ − 2
𝑑𝑡 = ∫ 2−𝑝
𝑑𝑡
𝑎
𝑥 1/𝑎
𝑡 0
𝑡
∞ 1/𝑎
1
∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑡−2 𝑓 ( ) 𝑑𝑡
𝑎 0
𝑡