Numerical Methods in Engineering (ENGR20005)
Lecture 11 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L11.1:
Numerical Integration: Trapezoidal Rule
2
“Book” (Section 6.2, pg. 109)
Integration
Integration means finding the area under a curve between two points, x0 and x1.
f(x1) f (x0)
x0
x1
3
Integration – Trapezoidal Rule
f (x1) f (x0)
f (x1) f (x0)
x0 x1 x0 x1
If we only have information at these two points, the best we can do is fit straight line through
these two points and find the area under it. We estimate A with A′ (see figure above).
A
A′
4
error
The error is quite large!
5
f (x1) f (x0)
Integration – Trapezoidal Rule
We estimate A with A′ (see figure below).
f (x1) f (x0)
x0 x1 x0 x1
We know that Lagrangian interpolation for two points can be written
as
A
A′
f(x)= (x−x1) f(x)+ (x−x0) f(x) 1 (x0 − x1) 0 (x1 − x0) 1
6
Integration – Trapezoidal Rule
f (x) = (x − x1) f (x ) + (x − x0) f (x ) 1 (x0 − x1) 0 (x1 − x0) 1
f (x1) f(x0)
A′
x0
To obtain area A′ all we have to do is to integrate f1(x) from x0 to x1.
∫x1 x0
f(x) ∫x1
0 (x − x1) d x
f1(x)d x =
+ f(x1) ∫x1(x−x0)dx
(x0−x1) x0 (x1−x0) x0
(11.1)
x1
7
Integration – Trapezoidal Rule
∫x1 f(x)dx= f(x0) ∫x1(x−x)dx x1 (x0−x1)x 1
(11.1)
0
f(x) ∫x0
1 1 (x − x )dx
+
∫ x1 f1(x)dx = x1 − x0 (f(x0) + f(x1))
(x1−x0) x 0 0
Integrate Eq. (11.1) to get
x0 2
∫ x1 f1(x)dx = Δ2 (f(x0) + f(x1))
(11.2)
x0
where Δ = x1 − x0. Eq. (11.2) is called the Trapezoidal rule:
In general, for any two points, xi and xi+1, the integral of any function can be approximated by the Trapezoidal rule as:
(11.3)
∫ xi+1 f(x)dx ≈ ∫ xi+1 f1(x)dx = Δ2 (f(xi) + f(xi+1)) xi xi
8
Integration – Trapezoidal Rule Error analysis
∫ xi+1 f(x)dx ≈ ∫ xi+1 f1(x)dx = Δ2 (f(xi) + f(xi+1)) xi xi
What is the error associated with the Trapezoidal rule?
∫ xi+1 xi
ε=− Exact integral, I.
Error
∫ xi+1 xi
f(x)dx
f(x)dx ≈ Δ2 (f(xi) + f(xi+1))
Δ2
(f(xi) + f(xi+1))
(11.4)
Approximation of I, using the trapezoidal rule.
9
ε=−
Perform Taylor series expansion of f(x) about x = xi
= =
where Δ = xi+1 − xi
Δ2 (f(xi) + f(xi+1)) =
Substitute Eqs. (11.6) and (11.5) into Eq. (11.4) to get, error
∫xi+1 Δ Δ3 Δ4
ε= x f(x)dx−2(f(xi)+f(xi+1))=−12f(2)(xi)−24f(3)(xi)+….
i
Therefore the local error is proportional to Δ3
10
(11.4)
∫xi+1 f(xi)+(x−xi)f′(xi)+ 1 (x−xi)2f(2)(xi)+ 1 (x−xi)3f(3)(xi)+…dx xi 2!3!
∫xi+1 f(x)dx xi
Δf(xi)+ Δ2 f′(xi)+ Δ3 f(2)(xi)+ Δ4 f(3)(xi)+…. 2! 3! 4!
(11.5)
Δ f(xi)+f(xi)+Δf′(xi)+ Δ2 f(2)(xi)+ Δ3 f(3)(xi)+… 2( 2! 3! )
(11.6)
∫ xi+1 xi
f(x)dx
Δ2
(f(xi) + f(xi+1))
Example L11.1: Use theTrapezoidal rule to approximate the following integrals.
I1 = ∫ 10 1x d x 2
I 2 = ∫ 26 1x d x
I3 = ∫24 1x d x
I 4 = ∫ 23 1x d x
Compare with the exact solution and show that the error gets smaller for smaller values of Δ
a) b ) c)
d )
11
I1 = ∫ 10 1x d x 2
2 10
I 2 = ∫ 26 1x d x
26
I 3 = ∫ 24 1x d x
I 4 = ∫ 23 1x d x
Smaller
error as Δ gets smaller
Integral Exact Trapezzoidal D abs(Error)
I1
1.60943791243 2.4 8 0.790562087565899
I2
1.09861228866 1.333333333 4 0.234721044665224
I3
0.6931471805 0.75 2 0.056852819440055
I4
0.4054651081 0.4166666666 1 0.011201558558502
24
23
12
L11.2:
Trapezoidal Rule: Large domain
13
“Book” (Section 6.2, pg. 109)
Integration – Trapezoidal Rule
≈
∫ xi+1 f1(x)dx = Δ2 (f(xi) + f(xi+1)) xi
∫xi+1 f(x)dx xi
f (x1) f (x0)
f (x1) f (x0)
A
A′
x0 Δx1
If size of interval, Δ, is large, the error is large as well.
x0 Δx1
14
Integration – Trapezoidal Rule
f (x1) f (x0)
error
f (x1) f (x0)
A
A′
Δ x1
x0 Δ x1
x0
If size of interval, Δ, is large, the error is large as well.
15
f (x1) f (x0)
Integration – Trapezoidal Rule
f (x1) f (x0)
A
A′
x0 Δx1 How can we minimize error as we integrate over a large interval?
Divide domain into smaller intervals, usually of equal length, apply the Trapezoidal rule in each of the subintervals
x0 Δx1
If size of interval, Δ, is large, the error is large as well.
16
Integration – Trapezoidal Rule
x0 Δ=Δ1 x1 x0 Δ=Δ1/2 x1 Δ=Δ1/2 x2 Intuitively, you get more accurate answers as you divide the large
interval up into smaller intervals. See figure above.
17
Integration – Trapezoidal Rule
‘negative’ error
x0 Δ=Δ1 x1 x0 Δ=Δ1/2 x1 Δ=Δ1/2 x2 Intuitively, you get more accurate answers as you divide the large
interval up into smaller intervals. See figure above.
18
Integration – Trapezoidal Rule
Less accurate approximation
More accurate approximation
x0 Δ=Δ1/2 x1 Δ=Δ1/2 x2
Intuitively, you get more accurate answers as you divide the large
interval up into smaller intervals. See figure above.
x0 x1 x2 x3 x4 Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4
19
Integration – Trapezoidal Rule
A1
A0 A2A3
x0 x1 x2 x3 x4
Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4
Intuitively, you get more accurate answers as you divide the large interval up into smaller intervals. See figure above.
What is general formula for Trapezoidal rule if we apply it to smaller subintervals ?
20
Integration – Trapezoidal Rule
Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4 Δ=Δ1/4 = Δ2 (f(x0) + f(x1)) + Δ2 (f(x1) + f(x2)) + Δ2 (f(x2) + f(x3)) + Δ2 (f(x3) + f(x4))
A = A0 + A1 + A2 + A3
= Δ2 (f(x0) + 2f(x1) + 2f(x2) + 2f(x3) + f(x4))
A0 A2A3
x0 x1 x2 x3 x4
= Δ2 (f(x0) + 2∑3 f(xi) + f(x4)) i=1
General formula for multiple application of trapezoidal rule to interval a ≤ x ≤ b, divided into N subintervals of equal size Δ, is
∫b a
f(x)dx ≈ Δ2 (f(a) + 2 21
N∑−1 i=1
f(xi) + f(b))
(11.7)
A1
We will now look at the error associated with the multiple application of Trapezoidal rule. Suppose we want to calculate the following integral between a ≤ x ≤ b.
∫ab f(x)dx The error associated with every interval is
εi=−Δ3 f(2)(xi)−Δ4 f(3)(xi)+…. 12 24
Hence, if you have N segments, the total error would be ∑N εi ≈∑N −Δ3 f(2)(xi)−Δ4 f(3)(xi)+….
i=1 i=1 12 24
If we ignore terms of the order higher than Δ4 and assume that all the f (2)(xi) are the same,
f (2)(xi) = f (2)(ξ) then
∑N ε i ≈ − N Δ 3 f ( 2 ) ( x i )
i=1 12 =NΔ3 f(2)(ξ)
RecallthatΔ= b−a N
(11.8)
12 3
=(b−a)Δ f(2)(ξ) Δ 12
=(b−a)Δ2 f(2)(ξ) 12
So the multiple-application of Trapezoidal rule has the error of the order of Δ2 22
Example L11.2: Use theTrapezoidal rule to approximate the following integral.
I1 = ∫ 10 1x d x 2
Use increasing numbers of intervals and compare with exact value of 1.6094379.
23
I1 = ∫ 10 1x d x 2
2 10
Exact Trapezzoidal D abs(Error)
1.60943791243 2.4 8 0.790562087565899
1.60943791243 1.866666666666 4 0.257228754232566
1.60943791243 1.683333333333 2 0.073895420899233
1.60943791243 1.6289682539 1 0.019530341534153
2 6 10
2 4 6 8 10
2 3 4 5 6 7 8 9 10
24
Lecture11M02.m
Delta=1 x=2:Delta:10 y=1./x;
I=MyTrapezoidal(y,Delta)
%Check with MATLAB Imatlab=Delta*trapz(y)
function I=MyTrapezoidal(y,dx) n=length(y);
end
∫b N−1
I=dx*(sum(y)-(y(1)+y(n))/2.0);
a
f(x)dx ≈ Δ2 (f(a) + 2 ∑ f(xi) + f(b)) i=1
25
L11.3: Simpson’s Rule
26
“Book” (Section 6.3, pg. 113)
f (x2) f (x0)
Integration – Simpson’s Rule
f (x2) f (x0)
f(x1)
x2 x0 x1 x2
x0
Trapezoidal rule: fit Lagrange polynomial through two points.
Simpson rule: Increase accuracy by fitting Lagrange polynomial through three equally spaced points and find area under Lagrange polynomial.
27
Integration – Simpson’s Rule
Error
The Lagrange interpolating Polynomial can be written as
f(x)= (x−x1)(x−x2) f(x)
f(x2) f(x0)
2
(x0 −x1)(x0 −x2) 0 + (x−x0)(x−x2) f(x )
(x1 − x0) (x1 − x2) 1 +(x−x0)(x−x1) f(x)
(x2 −x0)(x2 −x1) 2
f (x1)
x0
Δ x2
with x1 − x0 = x2 − x1 = Δ, integrate between x0 and x2:
Δ x1
x∫2 Δ
f2(x) = 3 (f(x0) + 4f(x1) + f(x2))
In general, for any three points, xi, xi+1 and xi+2, the integral of any function can be
x0 approximated by the Simpson rule as:
x∫
xi 28
i+2 Δ
f2(x) = 3 (f(xi) + 4f(xi+1) + f(xi+2))
Integration – Simpson’s Rule
In general, for any three points, xi, xi+1 and xi+2, the integral of any function can be approximated by the Simpson rule as:
x∫ xi
i+2 Δ( )
f2(x) = 3 f(xi) + 4f(xi+1) + f(xi+2) (Can be) More accurate than Trapezoidal rule.
(11.9)
29
Integration – Simpson Rule Error analysis
What is the error associated with the Simpson rule?
∫ xi+2 xi
ε=− Exact integral, I.
Error
∫ xi+2 xi
f(x)dx
f(x)dx ≈ Δ3 (f(xi) + 4f(xi+1) + f(xi+2))
Δ3
(f(xi) + 4f(xi+1) + f(xi+2))
Approximation of I, using the trapezoidal rule.
30
∫ xi+2 xi
f(x)dx
ε = − Δ3 (f(xi) + 4f(xi+1) + f(xi+2)) Do a Taylor series expansion of f(x) about x = xi
(11.10)
∫xi+2 f(x)dx xi
=∫xi+2 f(xi)+(x−xi)f′(xi)+ 1 (x−xi)2f(2)(xi)+ 1 (x−xi)3f(3)(xi)+…dx xi 2!3!
=
(11.11)
2Δf(xi)+ 2Δ2 f′(xi)+ 8Δ3 f(2)(xi)+ 16Δ4 f(3)(xi)+ 32Δ5 f(4)(xi)+…. 2! 3! 4! 5!
Δ3 (f(xi) + 4f(xi+1) + f(xi+2))
= Δ3 f ( x i ) +
4Δ f(xi)+Δf′(xi)+ Δ2 f(2)(xi)+ Δ3 f(3)(xi)+ Δ4 f(4)(xi)+…
3( 2! 3! 4! )
Δ f(xi)+2Δf′(xi)+ 4Δ f(2)(xi)+ 8Δ f(3)(xi)+ 16Δ f(4)(xi)+… 3( 2!2 3!3 4!4 )
Substitute Eqs (11.12) and (11.11) into Eq. (11.10) gives
∫xi+1 Δ Δ5
ε= x f(x)dx−3(f(xi)+4f(xi+1)+f(xi+2))=−90f(4)(xi)+….
i
Therefore the local error is proportional to Δ5
31
(11.12)
Example L11.3: Use the Simpson rule to approximate the following integrals.
I1 = ∫ 10 1x d x 2
I 2 = ∫ 26 1x d x
I3 = ∫24 1x d x
I 4 = ∫ 23 1x d x
Compare with the exact solution and show that the error gets smaller for smaller values of Δ
a) b ) c)
d )
32
Δ=4
I1 = ∫ 10 1x d x 2
10
I2 = ∫26 1x d x I 3 = ∫ 24 1x d x
Integral Exact Simpson D abs(Error)
I1
1.60943791243 1.6888888889 4 0.079450976454788
I2
1.09861228866 1.11111111111 2 0.012498822443002
I3
0.6931471805 0.6944444444 1 0.001297263884499
I4
0.4054651081 0.405555555 0.5 0.00009044744739
2
Smaller
error as Δ gets smaller
Δ = 2
26
24
I 4 = ∫ 23 1x d x 23
33
Multiple application of Simpson’s Rule
f (x)
∫x8 f(x)dx≈ Δ3 (f(x0)+4f(x1)+f(x2))+ Δ3 (f(x2)+4f(x3)+f(x4)) x0
+3 f(x)+4f(x)+f(x) +3 f(x)+4f(x)+f(x) Δ( 4 5 6) Δ( 6 7 8)
= Δ3[f(x0)+4(f(x1)+f(x3)+f(x5)+f(x7)) +2(f(x2)+f(x4)+f(x6))+f(x8)]
x0 x1 x2 x3 x4 x5 x6 x7 x8
In general, if you have N intervals, Simpson’s rule can be written as
∫x N−1 N−2
N f(x)dx ≈ Δ3 (f(x0)+4 ∑ f(xi)+2 ∑ f(xi)+f(xN))
x0 i=1iodd i=1ieven
(11.13)
Note that Eq. (11.13) only works if N is even i.e. the total number of data points (N + 1) is odd.
34
35
36
37
38
L11.4:
Romberg Integration: Richardson’s extrapolation
39
Richardson’s Extrapolation
The exact value of an integral can be decomposed into its numerical approximation and an error term.
I = I(Δ) + E(Δ)
(11.14)
I is the exact value of the integral
I(Δ) is the integral predicted by the integration scheme for spacing Δ E(Δ) is the error.
If we make two predictions of I, one with spacing Δ1 and also with spacing Δ2, then
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
We will assume that Δ2 < Δ1. If we look at the multiple application of the
Trapezoidal rule, we know from error analysis (see Eq. (11.8)) that
E(Δ1) ≈ Δ21 E(Δ2) Δ2
if we assume that the double derivative, f (2)(ξ) term is constant. 40
(11.15)
(11.16)
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
Substitute Eq. (11.16) into Eq. (11.15) gives
(11.15) E(Δ1) ≈ Δ21 E(Δ2) Δ2
(11.16)
I(Δ1)+E(Δ2)Δ21 ≈I(Δ2)+E(Δ2) Δ2
E(Δ2) ≈ I(Δ1) − I(Δ2) 1 − ( Δ1 )2
Substitute back into Eq. (11.15) gives
I = I(Δ2) + E(Δ2) ≈ I(Δ2) + I(Δ1) − I(Δ2) 1 − ( Δ1 )2
Δ2
Δ2
1
approximated values of the integral, I(Δ1) and I(Δ2). 41
−I(Δ1)
Trapezoidal rule. It gives you a much better estimate of the integral if you have two
I≈I(Δ2) 1+
Equation (11.17) is the Richardson extrapolation formula for the multiple application of the
( Δ1 )2 − 1 Δ2
( Δ1 )2 − 1 Δ2
(11.17)
1
Example L11.4: Use the Romberg Integration and theTrapezoidal rule to calculate integral:
I1 = ∫ 10 1x d x 2
Use Δ1 = 1.0 and Δ2 = 0.5 and compare with the exact value of 1.6094379.
I(Δ = 1.0)=1.628968253968254 I(Δ = 0.5)=1.614406323810349
From Eq. (11.17)
1+1
I ≈ I(Δ2) − I(Δ1) Δ1 =2
( Δ1 )2 − 1 Δ2
4 3
1
( Δ1 )2 − 1 Δ2
1 3
Δ2
I ≈ I(Δ = 0.5) 43 − I(Δ = 1.0) 13 = 1.609552347091047
42