Chapter 4 Numerical Integration Numerical integration is a primary tool used by engineers and scientists to obtain approximate answers for definite integrals that cannot be solved analytically. In the area of statistical thermodynamics, the Debye model for calculating the heat capacity of a solid involves the following function 更(x) Since there is no analytic expression for D(r), numerical integration must be used to obtain approximate values. For example, the value d(5)is the area under the curve Figure 4.1 The area under the curve y=f(t)for 0<t<5
Chapter 4 Numerical Integration Numerical integration is a primary tool used by engineers and scientists to obtain approximate answers for definite integrals that cannot be solved analytically. In the area of statistical thermodynamics, the Debye model for calculating the heat capacity of a solid involves the following function; Φ(x) = Z x 0 t 3 e t − 1 dt. Since there is no analytic expression for Φ(x), numerical integration must be used to obtain approximate values. For example, the value Φ(5) is the area under the curve 0 1 2 3 4 5 6 7 0 0.5 1 1.5 Figure 2.1 y=f(t) Figure 4.1 The area under the curve y = f(t) for 0 ≤ t ≤ 5. 4
Figure 4.1 Values of (a) 1.00.2248052 2.01.1763426 3025522185 4.03.8770542 504899892 6.055858554 7.06.0031690 8.06.2396238 9.06.3665739 1006.4319219 y=f(t=t/(e-1)for 0<t<5(see Figure 4.1). The numerical approximation for 9()2=a=1≈48892 Each additional value of (a)must be determined by another numerical integration Table 4. 1 lists several of these approximations over the interval [1,10 ti. The purpose of this chapter is to develop the basic principles of numerical integra- tion. In Chapter 9, numerical integration formulas are used to derive the predictor corrector methods for solving differential equations 4.1 Introduction to Quadrature We now approach the subject of numerical integration. The goal is to approximate the definite integral of f(a) over the interval a, b by evaluating f()at a finite number of ample point Definition 4.1 Suppose that a= to <1<...<M=b. A formula of the form QU]=∑mkf(xk)=tof(xo)+mf(x1)+…+Mf(xM) with the property that f(a d x=QU]+Elf (4.2)
Figure 4.1 Values of Φ(x) x Φ(x) 1.0 0.2248052 2.0 1.1763426 3.0 2.5522185 4.0 3.8770542 5.0 4.8998922 6.0 5.5858554 7.0 6.0031690 8.0 6.2396238 9.0 6.3665739 10.0 6.4319219 y = f(t) = t 3/(e t − 1) for 0 ≤ t ≤ 5 (see Figure 4.1). The numerical approximation for Φ(5) is Φ(5) = Z 5 0 t 3 e t − 1 dt ≈ 4.8998922. Each additional value of Φ(x) must be determined by another numerical integration. Table 4.1 lists several of these approximations over the interval [1, 10]. The purpose of this chapter is to develop the basic principles of numerical integration. In Chapter 9, numerical integration formulas are used to derive the predictorcorrector methods for solving differential equations. 4.1 Introduction to Quadrature We now approach the subject of numerical integration. The goal is to approximate the definite integral of f(x) over the interval [a, b] by evaluating f(x) at a finite number of sample points. Definition 4.1 Suppose that a = x0 < x1 < · · · < xM = b. A formula of the form Q[f] = X M k=0 wkf(xk) = w0f(x0) + w1f(x1) + · · · + wMf(xM) (4.1) with the property that Z b a f(x)dx = Q[f] + E[f] (4.2) 5
is called a numerical integration or quadrature formula. The term E() is called the truncation error for integration. The values ak k-o are called the quadrature nodes, and wk ko are called the weights taa epending on the application, the nodes ck) are chosen in various ways. For the pezoidal rule, Simpsons rule, and Boole's rule, the nodes are chosen to be equally spaced. For Gauss-Legendre quadrature, the nodes are chosen to be zeros of certain Legendre polynomials. When the integration formula is used to develop a predictor formula for differential equations, all the nodes are chosen less than b. For all applica- tion, it is necessary to know something about the accuracy of the numerical solution Definition 4.2. The degree of precision of a quadrature formula is the positive integer n such that E[P]=0 for all polynomials Pi(a) of degree i<n, but for which E[P+1 +0 for some polynomial Pn+1()of degree n+1 The form of E[Pl can be anticipated by studying what happens when f(a) polynomial. Consider the arbitrary polynomial ying P(x)=a1x2+x1-1x2-1+ +al.+ ao of degree i If< n, then Pin+()=0 for all x, and P(n+ ()=(n+1)lan-1 for all Thus it is not surprising that the general form for the truncation error term is Eff]=Kf(n+(c) where K is a suitably chosen constant and n is the degree of precision. The proof of this general result can be found in advanced books on numerical integration The derivation of quadrature formulas is sometimes based on polynomial interpo- lation. Recall that there exists a unique polynomial PM(a) of degree M passing through the M+l equally spaced points i(ck, yk)lMo. When this polynomial is used to approximate f(a), over a, bl, and then the integral of f(a)is approximated by the integral of PM(), the resulting formula is called a Newton-Cotes quadrature formula(see Figure 2.2). When the sample points o= a and M= b are used, it alled a closed Newton-Cotes formula. The next result gives the formulas when approximating polynomials of degree M= 1, 2, 3, and 4 are used Theorem 4.1(Closed Newton-Cotes Quadrature Formula). Assume that k o+ kh are equally spaced nodes and f&=f(Ck). The first four closed Newton-Cotes
is called a numerical integration or quadrature formula. The term E(f) is called the truncation error for integration. The values {xk} M k=0 are called the quadrature nodes, and {wk} M k=0 are called the weights. Depending on the application, the nodes {xk} are chosen in various ways. For the trapezoidal rule, Simpson’s rule, and Boole’s rule, the nodes are chosen to be equally spaced. For Gauss-Legendre quadrature, the nodes are chosen to be zeros of certain Legendre polynomials. When the integration formula is used to develop a predictor formula for differential equations, all the nodes are chosen less than b. For all application, it is necessary to know something about the accuracy of the numerical solution. Definition 4.2. The degree of precision of a quadrature formula is the positive integer n such that E[Pi ] = 0 for all polynomials Pi(x) of degree i ≤ n, but for which E[Pn+1] 6= 0 for some polynomial Pn+1(x) of degree n + 1. The form of E[Pi ] can be anticipated by studying what happens when f(x) is a polynomial. Consider the arbitrary polynomial Pi(x) = aix i + xi−1x i−1 + · · · + a1x + a0 of degree i. If i ≤ n, then P (n+1) i (x) ≡ 0 for all x, and P (n+1) n+1 (x) = (n + 1)!an−1 for all x. Thus it is not surprising that the general form for the truncation error term is E[f] = Kf(n+1)(c), (4.3) where K is a suitably chosen constant and n is the degree of precision. The proof of this general result can be found in advanced books on numerical integration. The derivation of quadrature formulas is sometimes based on polynomial interpolation. Recall that there exists a unique polynomial PM(x) of degree ≤ M passing through the M + 1 equally spaced points {(xk, yk)} M k=0. When this polynomial is used to approximate f(x), over [a, b], and then the integral of f(x) is approximated by the integral of PM(x), the resulting formula is called a Newton-Cotes quadrature formula (see Figure 2.2). When the sample points x0 = a and xM = b are used,it is called a closed Newton-Cotes formula. The next result gives the formulas when approximating polynomials of degree M = 1, 2, 3, and 4 are used. Theorem 4.1 (Closed Newton-Cotes Quadrature Formula). Assume that xk = x0 + kh are equally spaced nodes and fk = f(xk). The first four closed Newton-Cotes 6
quadrature formulas are ()d≈=2(6+/)( (the trapezoidal r f(x)dx≈n(后+4f1+f2)( Simpson' s rule), f(a)dr s o(o+3f1+3f2+f3)(Simpsons rule), 层f()≈(7+32h+125+32)5+7)( Boole's rule 4. Figure 2.2(a) The trapezoidal rule integrates( b) Simpson's rule integrates(c)Simp son's rule integrates(d)Boole's rule integrates Corollary 4.1(Newton-Cotes Precision). Assume that f(a)is sufficiently dif- ferentiable; then EU for Newton-Cotes quadrature involves an appropriate higher derivative. The trapezoidal rule has degree of precision n=1. If f E Ca, bl,then h3r(2(c) f(adr=o+f1)-iofo
quadrature formulas are Z x1 x0 f(x)dx ≈ h 2 (f0 + f1) (the traezoidal rule), (4.4) Z x2 x0 f(x)dx ≈ h 3 (f0 + 4f1 + f2) (Simpson’s rule), (4.5) Z x3 x0 f(x)dx ≈ 3h 8 (f0 + 3f1 + 3f2 + f3) (Simpson’s 3 8 rule), (4.6) Z x4 x0 f(x)dx ≈ 2h 45 (7f0 + 32f1 + 12f2 + 32f3 + 7f4) (Boole’s rule). (4.7) Figure 2.2 (a) The trapezoidal rule integrates (b) Simpson’s rule integrates (c) Simpson’s rule integrates (d) Boole’s rule integrates Corollary 4.1 (Newton-Cotes Precision). Assume that f(x) is sufficiently differentiable; then E[f] for Newton-Cotes quadrature involves an appropriate higher derivative. The trapezoidal rule has degree of precision n = 1. If f ∈ C 2 [a, b], then Z x1 x0 f(x)dx = h 2 (f0 + f1) − h 3 12 f (2)(c). (4.8) 7
Simpson's rule has degree of precision n=3. If f E C4a, b],then f(c)h f o +4f1 +f2)-oaf(c) (4.9) Simpson's& rule has degree of precision =3. If f E C [a,b],then ∫(x)d (f0+3f1+3f2+f3)-f4(c) Boole's rule has degree of precision n= 5, If f E C a, b],then 2h 8h ∫(x)d (7f0+32f1+12/2+32/3+7f4) 945(c) 4.11 Proof of Theorem 4. 1. Start with the Lagrange polynomial PM(ar)based on ro, T1,.,M that can be used to approximate f(a) f(x)≈PMr(x)=∑LMAk(x) (4.12) where f&=f(ak)for k=0, 1,..., M. An approximation for the integral is obtained by replacing the integrand f(a) with the polynomial PM(a). This is the general method for obtaining a Newton-Cotes integration formula ∫(x) fR.k(a)d (4.13) k(a)d rfK U The details for the general computations of coefficients of wk in(4.13)are tedious. We shall give a sample proof of Simpsons rule, which is the case M= 2. This case involves the approximating polynomial a-1(r-ar P2(x)=f0 2)c(x-xo)(x-x2),( +f1 x-x0)(x-x1) 2 (x-x0)(x-x2)12(x2-x0(2-)(414)
Simpson’s rule has degree of precision n = 3. If f ∈ C 4 [a, b], then Z x2 x0 f(x) = h 3 (f0 + 4f1 + f2) − h 5 90 f (4)(c). (4.9) Simpson’s 3 8 rule has degree of precisionn = 3. If f ∈ C 4 [a, b], then Z x3 x0 f(x)dx = 3h 8 (f0 + 3f1 + 3f2 + f3) − 3h 5 80 f (4)(c). (4.10) Boole’s rule has degree of precision n = 5, If f ∈ C 6 [a, b], then Z x4 x0 f(x)dx = 2h 45 (7f0 + 32f1 + 12f2 + 32f3 + 7f4) − 8h 7 945 f (6)(c). (4.11) Proof of Theorem 4.1. Start with the Lagrange polynomial PM(x) based on x0, x1, . . . , xM that can be used to approximate f(x): f(x) ≈ PM(x) = X M k=0 fkLM,k(x), (4.12) where fk = f(xk) for k = 0, 1, . . . , M. An approximation for the integral is obtained by replacing the integrand f(x) with the polynomial PM(x). This is the general method for obtaining a Newton-Cotes integration formula: Z xM x0 f(x) ≈ Z xM x0 PM(x)dx = Z xM x0 ÃX M k=0 fkLM,k(x) ! dx = X M k=0 µZ xM x0 fkLM,k(x)dx¶ (4.13) = X M k=0 µZ xM x0 LM,k(x)dx¶ fk = X M k=0 wkfk. The details for the general computations of coefficients of wk in (4.13) are tedious. We shall give a sample proof of Simpson’s rule, which is the case M = 2. This case involves the approximating polynomial P2(x) = f0 (x − x1)(x − x2) (x0 − x1)(x0 − x2) + f1 (x − x0)(x − x2) (x − x0)(x − x2) + f2 (x − x0)(x − x1) (x2 − x0)(x2 − x1) . (4.14) 8