Numerical Differentiation and Integration

師大

Numerical Differentiation and Integration

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

January 1, 2008

師大

Outline

1 Numerical Differentiation

2 Richardson Extrapolation Method

3 Elements of Numerical Integration Newton-Cotes Formulas

Composite Newton-Cotes Forumlas

師大

Numerical Differentiation

f⁰(x₀) = lim

h→0

f (x0+ h) − f (x0)

h .

Question How accurate is

f (x0+ h) − f (x0)

h ?

Suppose a given function f has continuous first derivative and f⁰⁰ exists.

From Taylor’s theorem

f (x + h) = f (x) + f⁰(x)h + 1

2f⁰⁰(ξ)h², where ξ is between x and x + h, one has

f⁰(x) = f (x + h) − f (x)

h − h

2f⁰⁰(ξ) = f (x + h) − f (x)

h + O(h).

師大

Hence it is reasonable to use the approximation f⁰(x) ≈ f (x + h) − f (x)

h

which is called forward finite difference, and the error involved is

|e| = h

2|f⁰⁰(ξ)| ≤ h 2 max

t∈(x,x+h)|f⁰⁰(t)|.

Similarly one can derive the backward finite difference approximation f⁰(x) ≈ f (x) − f (x − h)

h (1)

which has the same order of truncation error as the forward finite difference scheme.

師大

The forward difference is an O(h) scheme. An O(h²) scheme can also be derived from the Taylor’s theorem

f (x + h) = f (x) + f⁰(x)h +1

2f⁰⁰(x)h²+1

6f⁰⁰⁰(ξ₁)h³ f (x − h) = f (x) − f⁰(x)h +1

2f⁰⁰(x)h²−1

6f⁰⁰⁰(ξ₂)h³, where ξ1 is between x and x + h and ξ2 is between x and x − h. Hence

f (x + h) − f (x − h) = 2f⁰(x)h +1

6[f⁰⁰⁰(ξ₁) + f⁰⁰⁰(ξ₂)]h³ and

f⁰(x) = f (x + h) − f (x − h)

2h − 1

12[f⁰⁰⁰(ξ1) + f⁰⁰⁰(ξ2)]h² Let

M = max

z∈[x−h,x+h]f⁰⁰⁰(z) and m = min

z∈[x−h,x+h]f⁰⁰⁰(z).

師大

If f⁰⁰⁰ is continuous on [x − h, x + h], then by the intermediate value theorem, there exists ξ ∈ [x − h, x + h] such that

f⁰⁰⁰(ξ) = 1

2[f⁰⁰⁰(ξ1) + f⁰⁰⁰(ξ2)].

Hence

f⁰(x) = f (x + h) − f (x − h)

2h −1

6f⁰⁰⁰(ξ)h² = f (x + h) − f (x − h)

2h + O(h²).

This is called center difference approximation and the truncation error is

|e| = h² 6 f⁰⁰⁰(ξ)

Similarly, we can derive an O(h²) scheme from Taylor’s theorem for f⁰⁰(x) f⁰⁰(x) = f (x + h) − 2f (x) + f (x − h)

h² − 1

12f⁽⁴⁾(ξ)h², where ξ is between x − h and x + h.

師大

Polynomial Interpolation Method

Suppose that (x₀, f (x₀)), (x₁, f (x₁)) · · · , (x_n, f (x_n) have been given, we apply the Lagrange polynomial interpolation scheme to derive

P (x) =

n

X

i=0

f (x_i)L_i(x), where

Li(x) =

n

Y

j=0,j6=i

x − x_j xi− x_j. Since f (x) can be written as

f (x) =

n

X

i=0

f (xi)L(x) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξx)w(x), where

w(x) =

n

Y

j=0

(x − x_j),

師大

we have,

f⁰(x) = X

i=0

nf (xi)L⁰_i(x) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξx)w⁰(x)

+ 1

(n + 1)!w(x) d

dxf⁽ⁿ⁺¹⁾(ξx).

Note that

w⁰(x) =

n

X

j=0 n

Y

i=0,i6=j

(x − x_i).

Hence a reasonable approximation for the first derivative of f is f⁰(x) ≈

n

X

i=0

f (xi)L⁰_i(x).

When x = x_k for some 0 ≤ k ≤ n,

w(xk) = 0 and w⁰(xk) =

n

Y

i=0,i6=k

(xk− x_i).

師大

Hence f⁰(x_k) =

n

X

i=0

f (x_i)L⁰_i(x_k) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ_x)

n

Y

i=0,i6=k

(x_k− x_i), (2)

which is called an (n + 1)-point formula to approximate f⁰(x).

• Three Point Formulas Since

L₀(x) = (x − x₁)(x − x₂) (x0− x₁)(x0− x₂) we have

L⁰₀(x) = 2x − x₁− x₂ (x0− x₁)(x0− x₂). Similarly,

L⁰₁(x) = 2x − x0− x₂

(x₁− x₀)(x₁− x₂) and L⁰₂(x) = 2x − x0− x₁ (x₂− x₀)(x₂− x₁).

師大

Hence

f⁰(xj) = f (x0)

 2x_j − x₁− x₂ (x0− x₁)(x0− x₂)



+ f (x1)

 2x_j− x₀− x₂ (x1− x₀)(x1− x₂)



+ f (x2)

 2xj − x₀− x₁ (x₂− x₀)(x₂− x₁)

 +1

6f⁽³⁾(ξj)

2

Y

k=0,k6=j

(xj− x_k),

for each j = 0, 1, 2. Assume that

x1 = x0+ h and x2 = x0+ 2h, for some h 6= 0.

Then

f⁰(x₀) = 1 h



−3

2f (x₀) + 2f (x₁) −1 2f (x₂)

 +h²

3 f⁽³⁾(ξ₀), f⁰(x₁) = 1

h



−1

2f (x₀) +1 2f (x₂)



−h²

6 f⁽³⁾(ξ₁), f⁰(x2) = 1

h

 1

2f (x0) − 2f (x1) +3 2f (x2)

 +h²

3 f⁽³⁾(ξ2).

師大

That is

f⁰(x0) = 1 h



−3

2f (x0) + 2f (x0+ h) −1

2f (x0+ 2h)

 +h²

3 f⁽³⁾(ξ0), f⁰(x0+ h) = 1

h



−1

2f (x0) +1

2f (x0+ 2h)



−h²

6 f⁽³⁾(ξ1), (3) f⁰(x₀+ 2h) = 1

h

 1

2f (x₀) − 2f (x₀+ h) +3

2f (x₀+ 2h)

 +h²

3 f⁽³⁾(ξ₂).(4) Using the variable substitution x0 for x0+ h and x0+ 2h in (3) and (4), respectively, we have

f⁰(x₀) = 1

2h[−3f (x₀) + 4f (x₀+ h) − f (x₀+ 2h)] +h²

3 f⁽³⁾(ξ₀),(5) f⁰(x₀) = 1

2h[−f (x₀− h) + f (x₀+ h)] − h²

6 f⁽³⁾(ξ₁), f⁰(x0) = 1

2h[f (x0− 2h) − 4f (x₀− h) + 3f (x₀)] + h²

3 f⁽³⁾(ξ2). (6) Note that (6) can be obtained from (5) by replacing h with −h.

師大

• Five-point Formulas

f⁰(x0) = 1

12h[f (x0− 2h) − 8f (x₀− h) + 8f (x₀+ h) − f (x0+ 2h)]

+ h⁴

30f⁽⁵⁾(ξ),

where ξ ∈ (x0− 2h, x₀+ 2h) and f⁰(x0) = 1

12h[−25f (x0) + 48f (x0+ h) − 36f (x0+ 2h) + 16f (x₀+ 3h) − 3f (x₀+ 4h)] +h⁴

5 f⁽⁵⁾(ξ), where ξ ∈ (x0, x0+ 4h).

師大

Round-off Error

Consider

f⁰(x₀) = 1

2h[−f (x₀− h) + f (x₀+ h)] − h²

6 f⁽³⁾(ξ₁),

where ^h₆²f⁽³⁾(ξ₁) is called truncation error. Let ˜f (x₀+ h) and ˜f (x₀− h) be the computed values of f (x0+ h) and f (x0− h), respectively. Then

f (x0+ h) = ˜f (x0+ h) + e(x0+ h) and

f (x₀− h) = ˜f (x₀− h) + e(x₀− h).

Therefore, the total error in the approximation f⁰(x₀) −f (x˜ ₀+ h) − ˜f (x₀− h)

2h = e(x₀+ h) − e(x₀− h)

2h − h²

6 f⁽³⁾(ξ₁) is due in part to round-off error and in the part to truncation error.

師大

Assume that

|e(x₀± h)| ≤ ε and |f⁽³⁾(ξ1)| ≤ M.

Then     

f⁰(x0) −f (x˜ 0+ h) − ˜f (x0− h) 2h

≤ ε 2h+h²

6 M ≡ e(h).

Note that e(h) attains its minimum at h =p3ε/M .³

師大

Richardson’s Extrapolation

Suppose ∀h 6= 0 we have a formula N₁(h) that approximates an unknown value M

M − N1(h) = K1h + K2h²+ K3h³+ · · · , (7) for some unknown constants K₁, K₂, K₃, . . .. If K₁ 6= 0, then the

truncation error is O(h). For example, f⁰(x) − f (x + h) − f (x)

h = −f⁰⁰(x)

2! h −f⁰⁰⁰(x)

3! h²−f⁽⁴⁾(x)

4! h³− · · · . Goal

Find an easy way to produce formulas with a higher-order truncation error.

Replacing h in (7) by h/2, we have M = N1

 h 2

 + K1

h 2 + K2

h² 4 + K3

h³

8 + · · · . (8)

師大

Subtracting (7) with twice (8), we get M = N2(h) −K2

2 h²−3K3

4 h³− · · · , (9)

where

N₂(h) = 2N₁ h 2



− N₁(h) = N₁ h 2

 +

 N₁ h

2



− N₁(h)

 , which is an O(h²) approximation formula.

Replacing h in (9) by h/2, we get M = N₂ h

2



−k₂

8 h²−3K₃

32 h³− · · · . (10) Subtracting (9) from 4 times (10) gives

3M = N₂ h 2



− N₂(h) +3K₃

8 h³+ · · · , which implies that

M =

 N2

 h 2



+N2(h/2) − N2(h) 3

 +K3

8 h³+ · · · ≡ N3(h) +K3

8 h³+ · · · .

師大

Using induction, M can be approximated by M = Nm(h) + O(h^m), where

N_m(h) = N_m−1 h 2



+N_m−1(h/2) − N_m−1(h) 2^m−1− 1 . Centered difference formula. From the Taylor’s theorem f (x + h) = f (x) + hf⁰(x) +h²

2!f⁰⁰(x) + h³

3!f⁰⁰⁰(x) + h⁴

4!f⁽⁴⁾(x) + h⁵

5!f⁽⁵⁾(x) + h⁶

6!f⁽⁶⁾(x) + · · · f (x − h) = f (x) − hf⁰(x) +h²

2!f⁰⁰(x) − h³

3!f⁰⁰⁰(x) + h⁴

4!f⁽⁴⁾(x) − h⁵

5!f⁽⁵⁾(x) + h⁶

6!f⁽⁶⁾(x) + · · · we have

f (x + h) − f (x − h) = 2hf⁰(x) + 2h³

3! f⁰⁰⁰(x) +2h⁵

5! f⁽⁵⁾(x) + · · · ,

師大

and, consequently,

f⁰(x0) = f (x0+ h) − f (x0− h)

2h − h²

3!f⁰⁰⁰(x0) +h⁴

5!f⁽⁵⁾(x0) + · · ·

 ,

≡ N₁(h) − h²

3!f⁰⁰⁰(x₀) +h⁴

5!f⁽⁵⁾(x₀) + · · ·



. (11)

Replacing h in (11) by h/2 gives f⁰(x₀) = N₁ h

2



−h²

24f⁰⁰⁰(x₀) − h⁴

1920f⁽⁵⁾(x₀) − · · · . (12) Subtracting (11) from 4 times (12) gives

f⁰(x0) = N2(h) + h⁴

480f⁽⁵⁾(x0) + · · · , where

N₂(h) = 1 3



4N₁ h 2



− N₁(h)



= N₁ h 2



+ N₁(h/2) − N₁(h)

3 .

師大

In general,

f⁰(x₀) = N_j(h) + O(h^2j) with

N_j(h) = N_j−1 h 2



+N_j−1(h/2) − N_j−1(h) 4^j−1− 1 .

Example

Suppose that x₀ = 2.0, h = 0.2 and f (x) = xe^x. Compute an

approximated value of f⁰(2.0) = 22.16716829679195 to six decimal places.

Solution. By centered difference formula, we have N1(0.2) = f (0.2 + h) − f (0.2 − h)

2h = 22.414160,

N₁(0.1) = f (0.1 + h) − f (0.1 − h)

2h = 22.228786.

師大

It implies that

N₂(0.2) = N₁(0.1) +N₁(0.1) − N₁(0.2)

3 = 22.166995

which does not have six decimal digits. Adding N1(0.05) = 22.182564, we get

N2(0.1) = N1(0.05) +N1(0.05) − N1(0.1)

3 = 22.167157

and

N₃(0.2) = N₂(0.1) +N₂(0.1) − N₂(0.2)

15 = 22.167168

which contains six decimal digits.

師大

O(h) O(h²) O(h³) O(h⁴) 1: N₁(h) = N (h)

2: N1(h/2) = N (h/2) 3: N2(h)

4: N₁(h/4) = N (h/4) 5: N₂(h/2) 6: N₃(h)

7: N₁(h/8) = N (h/8) 8: N₂(h/4) 9: N₃(h/2) 10: N₄(h)

師大

Elements of Numerical Integration

The basic method involved in approximating the integration Z b

a

f (x) dx, (13)

is called numerical quadrature and uses a sum of the type Z b

a

f (x) dx ≈

n

X

i=0

c_if (x_i). (14)

The method of quadrature in this section is based on the polynomial interpolation. We first select a set of distinct nodes {x0, x1, . . . , xn} from the interval [a, b]. Then the Lagrange polynomial

P_n(x) =

n

X

i=0

f (x_i)L_i(x) =

n

X

i=0

f (x_i)

n

Y

j=0

j6=i

x − xj

x_i− x_j

is used to approximate f (x). With the error term we have

師大

f (x) = P_n(x) + E_n(x) =

n

X

i=0

f (x_i)L_i(x) + f⁽ⁿ⁺¹⁾(ζx) (n + 1)!

n

Y

i=0

(x − x_i), where ζ_x∈ [a, b] and depends on x, and

Z _b

a

f (x) dx = Z _b

a

P_n(x) dx + Z _b

a

E_n(x) dx

=

n

X

i=0

f (x_i) Z b

a

L_i(x) dx + 1 (n + 1)!

Z b a

f⁽ⁿ⁺¹⁾(ζ_x)

n

Y

i=0

(x − x_i) dx.(15) The quadrature formula is, therefore,

Z b a

f (x) dx ≈ Z b

a

Pn(x) dx =

n

X

i=0

f (xi) Z b

a

Li(x) dx ≡

n

X

i=0

cif (xi), (16) where

ci = Z b

a

Li(x) dx = Z b

a n

Y

j=0

j6=i

x − xj

x_i− x_j dx. (17)

師大

Moreover, the error in the quadrature formula is given by

E = 1

(n + 1)!

Z _b

a

f⁽ⁿ⁺¹⁾(ζ_x)

n

Y

i=0

(x − x_i) dx, (18) for some ζ_x ∈ [a, b].

Let us consider formulas produced by using first and second Lagrange polynomials with equally spaced nodes. This gives the Trapezoidal rule and Simpson’s rule, respectively.

Trapezoidal rule: Let x₀ = a, x₁ = b, h = b − a and use the linear Lagrange polynomial:

P1(x) = (x − x1)

(x₀− x₁)f (x0) + (x − x0) (x₁− x₀)f (x1).

Then Z b

a

f (x)dx = Z x1

x0

 (x − x₁)

(x0− x₁)f (x0) + (x − x0) (x1− x₀)f (x1)

 dx +1

2 Z x1

x0

f⁰⁰(ζ(x))(x − x0)(x − x1)dx. (19)

師大

Theorem (Weighted Mean Value Theorem for Integrals) Suppose f ∈ C[a, b], the Riemann integral of g(x)

Z _b

a

g(x)dx = lim

max 4xi→0 n

X

i=1

g(x_i)4x_i,

exists and g(x) does not change sign on [a, b]. Then ∃ c ∈ (a, b) with Z b

a

f (x)g(x)dx = f (c) Z b

a

g(x)dx.

Since (x − x₀)(x − x₁) does not change sign on [x₀, x₁], by the Weighted Mean Value Theorem, ∃ ζ ∈ (x₀, x₁) such that

Z x1

x0

f⁰⁰(ζ(x))(x − x₀)(x − x₁)dx = f⁰⁰(ζ) Z x1

x0

(x − x₀)(x − x₁)dx

= f⁰⁰(ζ) x³

3 − x1+ x0

2 x²+ x0x1x

x1

x0

= −h³ 6 f⁰⁰(ζ).

師大

Consequently, Eq. (19) implies that Z b

a

f (x)dx =  (x − x₁)²

2(x₀− x₁)f (x0) + (x − x0)² 2(x₁− x₀)f (x1)

x1

x0

−h³ 12f⁰⁰(ζ)

= x₁− x₀

2 [f (x₀) + f (x₁)] − h³ 12f⁰⁰(ζ)

= h

2[f (x₀) + f (x₁)] − h³ 12f⁰⁰(ζ), which is called the Trapezoidal rule.

師大

If we choose x₀= a, x₁= ¹₂(a + b), x₂ = b, h = (b − a)/2, and the second order Lagrange polynomial

P₂(x) = f (x₀) (x − x₁)(x − x₂)

(x₀− x₁)(x₀− x₂) + f (x₁) (x − x₀)(x − x₂) (x₁− x₀)(x₁− x₂) +f (x2) (x − x0)(x − x1)

(x₂− x₀)(x₂− x₁) to interpolate f (x), then

師大

Z b a

f (x)dx = Z x2

x0

 (x − x1)(x − x2)

(x₀− x₁)(x₀− x₂)f (x0) + (x − x0)(x − x2) (x₁− x₀)(x₁− x₂)f (x1) + (x − x₀)(x − x₁)

(x₂− x₀)(x₂− x₁)f (x₂)

 dx +

Z _x2

x0

(x − x0)(x − x1)(x − x2)

6 f⁽³⁾(ζ(x))dx.

Since, letting x = x0+ th, Z x2

x0

(x − x₁)(x − x₂)

(x0− x₁)(x0− x₂)dx = h Z 2

0

t − 1

0 − 1 · t − 2 0 − 2dt

= h

2 Z 2

0

(t²− 3t + 2) dt = h 3, Z x2

x0

(x − x₀)(x − x₂)

(x1− x₀)(x1− x₂)dx = h Z 2

0

t − 0

1 − 0 · t − 2 1 − 2dt

= −h Z ₂

0

(t²− 2t) dt = 4h 3 ,

師大

Z x2

x0

(x − x₀)(x − x₁)

(x₂− x₀)(x₂− x₁)dx = h Z 2

0

t − 0

2 − 0 · t − 1 2 − 1dt

= h

2 Z 2

0

(t²− t) dt = h 3, it implies that

Z _b

a

f (x)dx = h 1

3f (x₀) +4

3f (x₁) +1 3f (x₂)



+ Z x2

x0

(x − x0)(x − x1)(x − x2)

6 f⁽³⁾(ζ(x))dx, which is called the Simpson’s rule.

Deriving Simpson’s rule in this way, however, provides only an O(h⁴) error term involving f⁽³⁾. A higher order error analysis can be derived by expanding f in the third Taylor’s formula about x1. Then for each x ∈ [a, b], there exists ζ_x∈ (a, b) such that

f (x) = f (x1) + f⁰(x1)(x − x1) + f⁰⁰(x1)

2 (x − x1)² +f⁰⁰⁰(x1)

6 (x − x1)³+f⁽⁴⁾(ζx)

24 (x − x1)⁴.

T.M. Huang (Nat. Taiwan Normal Univ.) Numerical Diff. & integ. January 1, 2008 29 / 31

師大

Then Z b

a

f (x) dx =



f (x1)(x − x1) +f⁰(x1)

2 (x − x1)²+f⁰⁰(x1)

6 (x − x1)³ +f⁰⁰⁰(x₁)

24 (x − x₁)⁴

    

b

a

+ 1 24

Z b a

f⁽⁴⁾(ζ_x)(x − x₁)⁴dx.

Note that (b − x₁) = h, (a − x₁) = −h, and since (x − x₁)⁴ does not change sign in [a, b], by the Weighted Mean-Value Theorem for Integral, there exists ξ1 ∈ (a, b) such that

Z b a

f⁽⁴⁾(ζx)(x − x1)⁴dx = f⁽⁴⁾(ξ1) Z b

a

(x − x1)⁴dx = 2f⁽⁴⁾(ξ1) 5 h⁵. Consequently,

Z b a

f (x) dx = 2f (x1)h +f⁰⁰(x1)

3 h³+f⁽⁴⁾(ξ1) 60 h⁵.

師大

Finally we replace f⁰⁰(x1) by the central finite difference formulation f⁰⁰(x1) = f (x0) − 2f (x1) + f (x2)

h² −f⁽⁴⁾(ξ2) 12 h², for some ξ₂ ∈ (a, b), to obtain

Z b a

f (x) dx = 2hf (x1) +h

3(f (x0) − 2f (x1) + f (x2))

−f⁽⁴⁾(ξ2)

36 h⁵+f⁽⁴⁾(ξ1) 60 h⁵

= h 1

3f (x0) +4

3f (x1) +1 3f (x2)



+ 1 90

 3

2f⁽⁴⁾(ξ1) −5

2f⁽⁴⁾(ξ2)

 h⁵. It can show that there exists ξ ∈ (a, b) such that

Z b a

f (x) dx = h

3 [f (x₀) + 4f (x₁) + f (x₂)] − f⁽⁴⁾(ξ) 90 h⁵. This gives the Simpson’s rule formulation.

師大

If |f⁽ⁿ⁺¹⁾(x)| ≤ M on [a, b], then 

Z b a

f (x) dx −

n

X

i=0

cif (xi)     

≤ M

(n + 1)!

Z b a

n

Y

i=0

(x − xi) dx. (20) The choice of nodes that makes the right-hand side of this error bound as small as possible is know to be

x_i = a + b

2 +b − a

2 cos (i + 1)π n + 2



, i = 0, 1, . . . , n. (21) Of course, a polynomial interpolation to f can be obtained in other ways, for example, polynomial in Newton’s form using divided-difference method,

Pn(x) = f (x0) +

n

X

i=1

f [x0, x1, . . . , xi]

i−1

Y

j=0

(x − xj)

where f [x0, x1, . . . , xi] are evaluated with the divided difference algorithm.

Then Z b

a

f (x) dx ≈ f (x₀)(b − a) +

n

Xf [x₀, x₁, . . . , x_i] Z b

a i−1

Y(x − x_j) dx. (22)

師大

The standard derivation of quadrature error formulas is based on determining the class of polynomials for which theses formulas produce exact results. The next definition is used to facilitate the discussion of this derivation.

Definition

The degree of accuracy, or precision, of a quadrature formula is the largest positive integer n such that the formula is exact for x^k, when

k = 0, 1, . . . , n.

The definition implies that the degree of accuracy of a quadrature formula is n if and only if the error E = 0 for all polynomials P (x) of degree less than or equal to n, but E 6= 0 for some polynomials of degree greater than n.

A quadrature formula of the form (14) is called a Newton-Cotes formula if the nodes {x0, x1, . . . , xn} are equally spaced. Consider a uniform

partition of the closed interval [a, b] by

x_i= a + ih, i = 0, 1, . . . , n, h = b − a n ,

師大

where n is a positive integer and h is called the step length. By introduction a new variable t such that x = a + ht, the fundamental Lagrange polynomial becomes

Li(x) =

n

Y

j=0

j6=i

x − xj

x_i− x_j =

n

Y

j=0

j6=i

a + ht − a − jh a + ih − a − jh =

n

Y

j=0

j6=i

t − j

i − j ≡ ϕ_i(t).

Therefore, the integration (17) gives

c_i = Z b

a

L_i(x) dx = Z n

0

ϕ_i(t)h dt = h Z n

0 n

Y

j=0

j6=i

t − j

i − j dt, (23)

and the general Newton-Cotes formula has the form Z b

a

f (x) dx = h

n

X

i=0

f (x_i) Z n

0 n

Y

j=0

j6=i

t − j

i − j dt+ 1 (n + 1)!

Z b a

f⁽ⁿ⁺¹⁾(ζ_x)

n

Y

i=0

(x−x_i) dx.

(24)

師大

The simplest case is to choose n = 1, x₀ = a, x₁ = b, h = b − a, and use the linear Lagrange polynomial

P₁(x) = f (x₀)x − x₁

x₀− x₁ + f (x₁) x − x₀

x₁− x₀ = f (a)x − b

a − b + f (b)x − a b − a. to interpolate f (x). Then

c0= h Z 1

0

t − 1

0 − 1dt = h

2, c1 = h Z 1

0

t − 0

1 − 0dt = h 2, and

Z b a

P₁(x) dx = c₀f (x₀) + c₁f (x₁) = h

2[f (a) + f (b)] .

Since (x − x₀)(x − x₁) = (x − a)(x − b) does not change sign on [a, b], by the Weighted Mean-Value Theorem for integrals, there exists some

師大

ξ ∈ (a, b) such that Z _b

a

f⁰⁰(ζ_x)(x − x₀)(x − x₁) dx = f⁰⁰(ξ) Z _b

a

(x − x₀)(x − x₁) dx

= f⁰⁰(ξ) Z b

a

(x − a)(x − b) dx

= f⁰⁰(ξ) 1 3x³−1

2(a + b)x²+ abx

    

b

a

= −1

6f⁰⁰(ξ)(b − a)³= −1

6f⁰⁰(ξ)h³. Consequently,

Z b a

f (x) dx = h

2 [f (a) + f (b)] −h³ 12f⁰⁰(ξ).

This gives the so-called Trapezoidal rule.

Trapezoidal Rule:

Z b

f (x) dx = 1

(b − a) [f (a) + f (b)] − h³

f⁰⁰(ξ), (25)

師大

where h = b − a and ξ ∈ (a, b).

It is evident that the error term of the Trapezoidal rule is O(h³). Since the rule involves f⁰⁰, it gives the exact result when applied to any function whose second derivative is identically zero, e.g., any polynomial of degree 1 or less. Hence the degree of accuracy of Trapezoidal rule is one.

If we choose n = 2, x0 = a, x1 = ¹₂(a + b), x2 = b, h = (b − a)/2, and the second order Lagrange polynomial

P₂(x) = f (x₀) (x − x₁)(x − x₂)

(x₀− x₁)(x₀− x₂)+f (x₁) (x − x₀)(x − x₂)

(x₁− x₀)(x₁− x₂)+f (x₂) (x − x₀)(x − x₁) (x₂− x₀)(x₂− x₁) to interpolate f (x), then

c0 = h Z 2

0

t − 1

0 − 1 · t − 2

0 − 2dt = h 2

Z 2 0

(t²− 3t + 2) dt = h 3, c1 = h

Z 2 0

t − 0

1 − 0 · t − 2

1 − 2dt = −h Z 2

0

(t²− 2t) dt = 4h 3 , c2 = h

Z 2 0

t − 0

2 − 0 · t − 1

2 − 1dt = h 2

Z 2 0

(t²− t) dt = h 3,

師大

and Z b

a

P2(x) dx = c0f (x0)+c1f (x1)+c2f (x2) = h 1

3f (a) +4

3f (a + b 2 ) +1

3f (b)



gives the so-called Simpson’s rule. Deriving the formulation this way, however, the error term

1 6

Z b a

f⁽³⁾(ζx)(x − x0)(x − x1)(x − x2) dx

provides only an O(h⁴) formulation involving f⁽³⁾. A higher order error analysis can be derived by expanding f in the third Taylor’s formula about x1. Then for each x ∈ [a, b], there exists ζx∈ (a, b) such that

f (x) = f (x₁)+f⁰(x₁)(x−x₁)+f⁰⁰(x₁)

(x−x₁)²+f⁰⁰⁰(x₁)

(x−x₁)³+f⁽⁴⁾(ζ_x)

(x−x₁)⁴.

師大

Then Z b

a

f (x) dx =



f (x₁)(x − x₁) +f⁰(x₁)

2 (x − x₁)²+f⁰⁰(x₁)

6 (x − x₁)³+ f⁰⁰⁰(x₁)

24 (x − x₁)⁴

    

b

a

+ 1 24

Z b a

f⁽⁴⁾(ζ_x)(x − x₁)⁴dx.

Note that (b − x1) = h, (a − x1) = −h, and since (x − x1)⁴ does not change sign in [a, b], by the Weighted Mean-Value Theorem for Integral, there exists ξ₁ ∈ (a, b) such that

Z _b

a

f⁽⁴⁾(ζ_x)(x − x₁)⁴dx = f⁽⁴⁾(ξ₁) Z _b

a

(x − x₁)⁴dx = 2f⁽⁴⁾(ξ1) 5 h⁵. Consequently,

Z b a

f (x) dx = 2f (x1)h +f⁰⁰(x1)

3 h³+f⁽⁴⁾(ξ1) 60 h⁵. Finally we replace f⁰⁰(x1) by the central finite difference formulation

f⁰⁰(x₁) = f (x₀) − 2f (x₁) + f (x₂)

h² −f⁽⁴⁾(ξ₂) 12 h²,

師大

for some ξ₂ ∈ (a, b), to obtain Z b

a

f (x) dx = 2hf (x1) +h

3(f (x0) − 2f (x1) + f (x2)) − f⁽⁴⁾(ξ2)

36 h⁵+f⁽⁴⁾(ξ1) 60 h⁵

= h 1

3f (x₀) +4

3f (x₁) +1 3f (x₂)

 + 1

90

 3

2f⁽⁴⁾(ξ₁) −5

2f⁽⁴⁾(ξ₂)

 h⁵.

By letting f (x) = x⁴, one can show that there exists ξ ∈ (a, b) such that Z b

a

f (x) dx = h 1

3f (x₀) +4

3f (x₁) +1 3f (x₂)



+f⁽⁴⁾(ξ) 90 h⁵. This gives the Simpson’s rule formulation.

Simpson’s Rule:

Z b a

f (x) dx = b − a 2

  1

3f (a) +4

3f (a + b 2 ) +1

3f (b)



+f⁽⁴⁾(ξ)

90 h⁵, (26) for some ξ ∈ (a, b). The Simpson’s rule is an O(h⁵) scheme and the degree of accuracy is three.

師大

The Trapezoidal and Simpson’s rules are examples of a class of methods known as closed Newton-Cotes formula. The (n + 1)-point closed Newton-Cotes method uses nodes xi= a + ih, for i = 0, 1, . . . , n, where h = (b − a)/n. Note that both endpoints, a = x₀ and b = x_n, of the closed interval [a, b] are included as nodes. The following theorem details the Newton-Cotes formulas and the associated error analysis.

師大

Theorem (Closed Newton-Cotes Formulas)

For a given function f (x) and closed interval [a, b], the (n + 1)-point closed Newton-Cotes method uses nodes

x_i= a + ih, i = 0, 1, . . . , n, h = b − a n . If n is even and f ∈ Cⁿ⁺²[a, b], then

Z b a

f (x) dx = h

n

X

i=0

α_if (x_i) +hⁿ⁺³f⁽ⁿ⁺²⁾(ξ) (n + 2)!

Z n 0

t²(t − 1) · · · (t − n) dt, (27) and if n is odd and f ∈ Cⁿ⁺¹[a, b], then

Z b a

f (x) dx = h

n

X

i=0

αif (xi) +hⁿ⁺²f⁽ⁿ⁺¹⁾(ξ) (n + 1)!

Z n 0

t(t − 1) · · · (t − n) dt, (28) where ξ ∈ (a, b) and

αi = Z n

0 n

Y

j=0

j6=i

t − j

i − j dt, i = 0, 1, . . . , n. (29)

Consequently, the degree of accuracy is n + 1 when n is an even integer,

T.M. Huang (Nat. Taiwan Normal Univ.) Numerical Diff. & integ. January 1, 2008 31 / 31

師大

The weights α_i in the Newton-Cotes formula has the property

n

X

i=0

αi = n. (30)

This can be shown by applying the formula to f (x) = 1 with interpolating polynomial Pn(x) = 1. Let s be the common denominator of αi, that is,

α_i = σ_i

s (⇒ σ_i = sα_i)

such that σ_i are integers, then the formulation for approximating the definite integral can be expressed as

Z b a

f (x) dx ≈ h

n

X

i=0

α_if (x_i) = h s

n

X

i=0

σ_if (x_i). (31)

Some of the most common closed Newton-Cotes formulas with their error terms are listed in the following table.

師大

Name n s σi Error

Trapezoidal rule 1 2 1, 1 −₁₂¹ f⁽²⁾(ξ)h³

Simpson’s rule 2 3 1, 4, 1 −₉₀¹ f⁽⁴⁾(ξ)h⁵

3/8-rule 3 ⁸₃ 1, 3, 3, 1 −₈₀³ f⁽⁴⁾(ξ)h⁵

Milne’s rule 4 ⁴⁵₂ 7, 32, 12, 32, 7 −₉₄₅⁸ f⁽⁶⁾(ξ)h⁷ 5 ²⁸⁸₅ 19, 75, 50, 50, 75, 19 −₁₂₀₉₆²⁷⁵ f⁽⁶⁾(ξ)h⁷ Weddle’s rule 6 140 41, 216, 27, 272, 27, 216, 41 −₁₄₀₀⁹ f⁽⁸⁾(ξ)h⁹ Another class of Newton-Cotes formulas is the open Newton-Cotes

formulas in which the nodes

xi = x0+ ih, i = 0, 1, . . . , n, where x0 = a + h and h = b − a n + 2, are used. This implies that xn= b − h, and the endpoints, a and b, are not used. Hence we label a = x−1 and b = x_n+1. The following theorem summarizes the open Newton-Cotes formulas.

師大

Theorem (Open Newton-Cotes Formulas)

For a given function f (x) and closed interval [a, b], the (n + 1)-point open Newton-Cotes method uses nodes

x_i = x₀+ ih, i = 0, 1, . . . , n, where x₀ = a + h and h = b − a n + 2. If n is even and f ∈ Cⁿ⁺²[a, b], then

Z b a

f (x) dx = h

n

X

i=0

αif (xi) +hⁿ⁺³f⁽ⁿ⁺²⁾(ξ) (n + 2)!

Z n+1

−1

t²(t − 1) · · · (t − n) dt, (32) and if n is odd and f ∈ Cⁿ⁺¹[a, b], then

Z b a

f (x) dx = h

n

X

i=0

α_if (x_i) +hⁿ⁺²f⁽ⁿ⁺¹⁾(ξ) (n + 1)!

Z n+1

−1

t(t − 1) · · · (t − n) dt, (33) where ξ ∈ (a, b) and

αi = Z n+1

−1 n

Y

j=0

j6=i

t − j

i − j dt, i = 0, 1, . . . , n. (34)

Consequently, the degree of accuracy is n + 1 when n is an even integer, and n when n is an odd integer.

T.M. Huang (Nat. Taiwan Normal Univ.) Numerical Diff. & integ. January 1, 2008 31 / 31

師大

The simplest open Newton-Cotes formula is choosing n = 0 and only using the midpoint x₀ = ^a+b₂ . Then the coefficient and the error term can be computed easily as

α0= Z −1

1

dt = 2, and h³f⁰⁰(ξ) 2!

Z 1

−1

t²dt = 1

3f⁰⁰(ξ)h³.

These gives the so-called Midpoint rule or Rectangular rule.

Midpoint Rule:

Z _b

a

f (x) dx = 2hf (x₀) +1

3f⁰⁰(ξ)h³ = (b − a)f (a + b 2 ) +1

3f⁰⁰(ξ)h³, (35)

for some ξ ∈ (a, b).

Analogous to the closed Newton-Cotes formulas, we list some of the commonly used open Newton-Cotes formulas in the following table.

師大

Name n s σ_i Error

Midpoint rule 0 1 2 ¹₃f⁽²⁾(ξ)h³ 1 2 3, 3 ³₄f⁽²⁾(ξ)h³ 2 3 8, −4, 8 ¹⁴₄₅f⁽⁴⁾(ξ)h⁵ 3 24 55, 5, 5, 55 ₁₄₄⁹⁵f⁽⁴⁾(ξ)h⁵

It is obvious that the Newton-Cotes formulas are generally not suitable for numerical integration over large interval. Higher degree formulas would be required, and the coefficients in these formulas are difficult to obtain. Also the Newton-Cotes formulas which are based on polynomial interpolation would be inaccurate over a large interval because of the oscillatory nature of high-degree polynomials. Now we discuss a piecewise approach, called composite rule, to numerical integration over large interval that uses the low-order Newton-Cotes formulas.

A composite rule is one obtained by applying an integration formula for a single interval to each subinterval of a partitioned interval. To illustrate the procedure, we choose an even integer n and partition the interval [a, b]

into n subintervals by nodes x₀ < x₁< · · · < x_n= b, and apply Simpson’s

師大

rule on each consecutive pair of subintervals. With

h = b − a

n and x_j = a + jh, j = 0, 1, . . . , n,

we have on each interval [x_2j−2, x_2j],

Z _x2j

x2j−2

f (x) dx = h

3 [f (x_2j−2) + 4f (x_2j−1) + f (x_2j)] − h⁵

90f⁽⁴⁾(ξ_j),

for some ξ ∈ (x , x ), provided that f ∈ C⁴[a, b]. The composite rule

師大

is obtained by summing up over the entire interval, that is, Z b

a

f (x) dx =

n/2

X

j=1

Z x2j

x2j−2

f (x) dx

=

n/2

X

j=1

 h

3(f (x2j−2) + 4f (x2j−1) + f (x2j)) − h⁵

90f⁽⁴⁾(ξj)



= h

3[f (x0) + 4f (x1) + f (x2) + f (x2) + 4f (x3) + f (x4) + f (x4) + 4f (x5) + · · · + f (xn−2) + 4f (xn−1) + f (xn)] − h⁵

90

n/2

X

j=1

f⁽⁴⁾(ξj)

= h

3[f (x0) + 4f (x1) + 2f (x2) + 4f (x3) + 2f (x4) + 4f (x5) + · · · + 2f (xn−2) + 4f (xn−1) + f (xn)] − h⁵

90

n/2

X

j=1

f⁽⁴⁾(ξj)

= h

3



f (x0) + 4

n/2

X

j=1

f (x2j−1) + 2

(n/2)−1

X

j=1

f (x2j) + f (xn)



− h⁵ 90

n/2

X

j=1

f⁽⁴⁾(ξj).

T.M. Huang (Nat. Taiwan Normal Univ.) Numerical Diff. & integ. January 1, 2008 31 / 31

師大

To estimate the error associated with approximation, since f ∈ C⁴[a, b], we have, by the Extreme Value Theorem,

min

x∈[a,b]f⁽⁴⁾(x) ≤ f⁽⁴⁾(ξ_j) ≤ max

x∈[a,b]f⁽⁴⁾(x), for each ξ_j ∈ (x_2j−2, x_2j). Hence

n 2 min

x∈[a,b]f⁽⁴⁾(x) ≤

n/2

X

j=1

f⁽⁴⁾(ξj) ≤ n 2 max

x∈[a,b]f⁽⁴⁾(x), and

x∈[a,b]min f⁽⁴⁾(x) ≤ 2 n

n/2

X

j=1

f⁽⁴⁾(ξj) ≤ max

x∈[a,b]f⁽⁴⁾(x).

By the Intermediate Value Theorem, there exists µ ∈ (a, b) such that

f⁽⁴⁾(µ) = 2 n

n/2

Xf⁽⁴⁾(ξ_j).

師大

Thus, by replacing n = (b − a)/h,

n/2

X

j=1

f⁽⁴⁾(ξj) = n

2f⁽⁴⁾(µ) = b − a

2h f⁽⁴⁾(µ).

Consequently, the composite Simpson’s rule is derived.

Composite Simpson’s Rule:

Z b a

f (x) dx = h 3



f (a) + 4

n/2

X

j=1

f (x_2j−1) + 2

(n/2)−1

X

j=1

f (x_2j) + f (b)



−b − a

180 f⁽⁴⁾(µ)h⁴, (36)

where n is an even integer, h = (b − a)/n, xj = a + jh, for j = 0, 1, . . . , n, and some µ ∈ (a, b).

The composite Midpoint rule can be derived in a similar way, except the midpoint rule is applied on each subinterval [x2j−1, x2j] instead. That is,

Z x2j

x2j−2

f (x) dx = 2hf (x_2j−1) +h³

3 f⁰⁰(ξ_j), j = 1, 2, . . . ,n 2.

師大

Note that n must again be even. Consequently, Z b

a

f (x) dx = 2h

n/2

X

j=1

f (x2j−1) +h³ 3

n/2

X

j=1

f⁰⁰(ξj).

The error term can be written as

n/2

X

j=1

f⁰⁰(ξj) = n

2f⁰⁰(µ) = b − a 2h f⁰⁰(µ),

for some µ ∈ (a, b). Therefore, the composite Midpoint rule has the following formulation.

Composite Midpoint Rule:

Z b a

f (x) dx = 2h

n/2

X

j=1

f (x2j−1) −b − a

6 f⁰⁰(µ)h², (37) where n is an even integer, h = (b − a)/n, x_j = a + jh, for

j = 0, 1, . . . , n, and some µ ∈ (a, b).

師大

To derive the composite Trapezoidal rule, we partition the interval [a, b] by n equally spaced nodes a = x₀ < x₁ < · · · < x_n= b, where n can be either odd or even. We then apply the trapezoidal rule on each subinterval

師大

[x_j−1, x_j] and sum them up to obtain Z b

a

f (x) dx =

n

X

j=1

Z xj

xj−1

f (x) dx

=

n

X

j=1

 h

2(f (xj−1) + f (xj)) −h³ 12f⁰⁰(ξj)



= h

2[f (x0) + f (x1) + f (x1) + f (x2) + · · · + f (xn−1) + f (xn)] − h³ 12

n

X

j=1

f⁰⁰(ξj)

= h

2[f (x₀) + 2f (x₁) + 2f (x₂) + · · · + 2f (x_n−1) + f (x_n)] − h³ 12

n

X

j=1

f⁰⁰(ξ_j)

= h

2



f (a) +

n−1

X

j=1

f (x_j) + f (b)



−h³ 12

n

X

j=1

f⁰⁰(ξ_j)

= h

2



f (a) +

n−1

Xf (x_j) + f (b)



−b − a

12 f⁰⁰(µ)h²,

師大

where each ξ_j ∈ (x_j−1, x_j) and µ ∈ (a, b).

Composite Trapezoidal Rule Z _b

a

f (x) dx = h 2



f (a) +

n−1

X

j=1

f (x_j) + f (b)



− b − a

12 f⁰⁰(µ)h², (38) where n is an integer, h = (b − a)/n, xj = a + jh, for j = 0, 1, . . . , n, and some µ ∈ (a, b).