Chapter 5: Numerical Integration and Differentiation

1 Chapter 5: Numerical Integration and DifferentiationPART I: Numerical IntegrationNewton-Cotes Integration FormulasThe idea of Newton-Cotes formulas is to replace a complicated function or tabu-lated data with an approximating function that is easy to baf(x)dx bafn(x)dxwherefn(x) =a0+a1x+a2x2+..+ Trapezoidal RuleUsing the first order Taylor series to approximatef(x),I= baf(x)dx baf1(x)dxwheref1(x) =f(a) +f(b) f(a)b a(x a)1 ThenI ba[f(a) +f(b) f(a)b a(x a)]dx= (b a)f(b) +f(a)2 The trapezoidal rule is equivalent to approximating the area of the trapezoidalFigure 1: Graphical depiction of the trapezoidal ruleunder the straight line connectingf(a)andf(b). An estimate for the local trun-2cation error of a single application of the trapezoidal rule can be obtained usingTaylor series asEt= 112f ( )(b a)3where is a value :Use the trapezoidal rule to numerically integratef(x) = + 25xfroma= 0tob= :f(a) =f(0) = , andf(b) =f(2) = (b a)f(b) +f(a)2= (2 0) + true solution is 20f(x)dx= ( + )|20= ( 2 + 22) 0 = (x)is a linear function, using the trapezoidal rule gets the exact :Use the trapezoidal rule to numerically integratef(x) = + 25x+ 3x23froma= 0tob= :f(0) = , andf(2) = (b a)f(b) +f(a)2= (2 0) + true solution is 20f(x)dx= ( + +x3)|20= ( 2 + 22+ 23) 0 = relative error is| t|= 100% = trapezoidal rule:Using smaller Integration interval can reduce the approximation error.

2 We candivide the Integration interval fromatobinto a number of segments and applythe trapezoidal rule to each segment. Divide(a, b)intonsegments of equalwidth. ThenI= baf(x)dx= x1x0f(x)dx+ x2x1f(x)dx+..+ xnxn 1f(x)dxwherea=x0< x1< .. < xn=b, andxi xi 1=h=b an, fori= 1,2, .. , the Trapezoidal rule for each integral yieldsI hf(x0) +f(x1)2+hf(x1) +f(x2)2+..+hf(xn 1) +f(xn)2=h2[f(x0) + 2n 1 i=1f(xi) +f(xn)]= (b a)f(x0) + 2 n 1i=1f(xi) +f(xn)2nThe approximation error using the multiple trapezoidal rule is a sum of the indi-vidual errors, ,Et= n i=1h312f ( i) = n i=1(b a)312n3f ( i)5 Letf = ni=1f ( i)n. Then the approximate error isEt= (b a)312n2f Example:Use the 2-segment trapezoidal rule to numerically integratef(x) = + 25x+ 3x2froma= 0tob= :n= 2,h= (a b)/n= (2 0)/2 = (0) = ,f(1) = , andf(2) = (b a)f(0) + 2f(1) +f(2)2n= 2 + 2 + relative error is| t|= 100% = s RulesAside from using the trapezoidal rule with finer segmentation, another way toimprove the estimation accuracy is to use higher order 2: Illustration of (a) Simpson s 1/3 rule, and (b) Simpson s 3/8 ruleSimpson s1/3rule:Given function values at 3 points as(x0, f(x0)),(x1, f(x1)), and(x2, f(x2)), we8can estimatef(x)using Lagrange polynomial interpolation.

3 ThenI= baf(x)dx x2x0[(x x1)(x x2)(x0 x1)(x0 x2)f(x0)+(x x0)(x x2)(x1 x0)(x1 x2)f(x1) +(x x0)(x x1)(x2 x0)(x2 x1)f(x2)]dxWhena=x0,b=x2,(a+b)/2 =x1, andh= (b a)/2,I h3[f(x0) + 4f(x1) +f(x2)] = (b a)f(x0) + 4f(x1) +f(x2)6It can be proved that single segment application of Simpson s1/3rule has atruncation error ofEt= 190h5f(4)( )where is s1/3rule yields exact results for third order polynomials even thoughit is derived from :Use Simpson s1/3rule to integratef(x) = + 25x+ 3x2+ 8x3froma= 0tob= :f(0) = ,f(1) = , andf(2) = (b a)f(0) + 4f(1) +f(2)6= 2 + 4 + exact integral is 20f(x)dx= ( + +x3+2x4)|20= ( 2+ 22+23+2 24) 0 = :Use Simpson s1/3rule to integratef(x) = + 25x+ 3x2+ 2x4froma= 0tob= :f(0) = ,f(1) = , andf(2) = (b a)f(0) + 4f(1) +f(2)6= 2 + 4 + exact integral is 20f(x)dx= ( + +x3+ )|20= ( 2+ 22+23+ 25) 0 = relative error is| t|= = Simpson s1/3ruleDividing the Integration interval intonsegments of equal width, we haveI= x2x0f(x)dx+ x4x2f(x)dx+.

4 + xnxn 2f(x)dxwherea=x0< x1< .. < xn=b, andxi xi 1=h= (b a)/n, fori= 1,2, .. , n. Substituting the Simpson s1/3rule for each integral yieldsI 2hf(x0) + 4f(x1) +f(x2)6+ 2hf(x2) + 4f(x3) +f(x4)6+ + 2hf(xn 2) + 4f(xn 1) +f(xn)6= (b a)f(x0) + 4 n 1i=1,3,5f(xi) + 2 n 2j=2,4,6f(xj) +f(xn)3nNote thatnhas to be :Use 4-segment Simpson s1/3rule to integratef(x) = + 25x+ 3x2+ 2x4froma= 0tob= :n= 4,h= (b a)/n= (x0) =f(0) = ,f(x1) =f( ) = ,f(x2) =f(1) = ,f(x3) =11f( ) = , andf(x4) =f(2) = (b a)f(0) + 4f( ) + 4f( ) + 2f(1) +f(2)3 4= 2 + 4 + 4 + 2 + exact integral is 20f(x)dx= ( + +x3+ )|20= ( 2+ 22+23+ 25) 0 = relative error is| t|= = s3/8ruleThis is to use a third-order Lagrange polynomial to fit to four points off(x)andyieldsI 3h8[f(x0) + 3f(x1) + 3f(x2) +f(x3)]whereh= (b a)/3. The approximation error using this rule isEt= 380h5f(4)( ) = (b a)56480f(4)( )where is of EquationsNewton-Cotes algorithms for equationsCompare the following two Pseudocodes for multiple applications of the trape-zoidal 1: Algorithm for multiple applications of the trapezoidal rulefunction Trapm(h,n,f)sum=f0for i=1:n-1sum=sum+2*fiendsum=sum+fnTrapm=h* sum/2 Pseudocode 2: Algorithm for multiple application of the trapezoidal rule whenfunctionf(x)is availablefunction TrapEq(n,a,b)h=(b-a)/nx=asum=f(x)for i=1:n-113x=x+hsum=sum+2*f(x)endsum=sum+f (b)TraEq=(b-a)*sum/(2*n)Pseudocode 1 can be used when only a limited number of points are given or thefunction is available.

5 Pseudocode 2 is for the case where the analytical functionis available. The difference between the two pseudocodes is that in Pseudocode2 neigher the independent nor the dependent variable values are passed into thefunction via its argument as in Pseudocode 1. When the analytical function isavailable, the function values are computed using calls to the function beinganalyzed,f(x).4 Romberg IntegrationRomberg Integration is one technique that can improve the results of numericalintegration using error-correction s extrapolationuses two estimates of an integral to compute a third,more accurate estimate and error associated with a multiple-application trapezoidal rule14can be represented asI=I(h) +E(h)whereIis the exact value of the integral,I(h)is the approximation from ann-segment application of the trapezoidal rule with step sizeh= (b a)/n, andE(h)is the truncation error.

6 If we make two separate estimates using step sizesofh1andh2and have exact values for the error, thenI(h1) +E(h1) =I(h2) +E(h2)(1)The error of the multiple-application trapezoidal rule can be represented approx-imately asE b a12h2 f whereh= (b a)/n. Assuming that f is constant regardless of step size, wehaveE(h1)E(h2) h21h22 Then we haveE(h1) E(h2)(h1h2)2which can be substituted into (1):I(h1) +E(h2)(h1h2)2=I(h2) +E(h2)15 ThenE(h2)can be solved asE(h2) =I(h2) I(h1)(h1/h2)2 1 This estimate can then be substituted intoI=I(h2) +E(h2)to yield an improved estimate of the integral:I I(h2) +I(h2) I(h1)(h1/h2)2 1It can be shown that the error of this estimate isO(h4). Thus, we have combinedtwo trapezoidal rule estimates ofO(h2)to yield a new estimate ofO(h4). For thespecial case where the interval ish2=h1/2, we haveI I(h2) +122 1[I(h2) I(h1)] =43I(h2) 13I(h1)With two improved integrals ofO(h4)on the basis of three trapezoidal rule esti-mates, we can combine them to yield an even better value withO(h6).

7 The Romberg Integration algorithmhas the general form asIj,k 4k 1Ij+1,k 1 Ij,k 14k 1 116whereIj+1,k 1andIj,k 1are the more and less accurate integrals, respectively,andIj,kis the improved integral. The indexksignifies the level of the integra-tion, wherek= 1corresponds to the original trapezoidal rule estimates,k= 2corresponds toO(h4),k= 3toO(h6), and so :f(x) = + 25x 200x2+ 675x3 900x4+ 400x5, find baf(x)dx,a= 0,b= :True solution:I= (x)dx= ( + 25x 200x2+ 675x3 900x4+400x5)dx= 1,h1= (0) = ,f( ) = ,1= (x)dx= f(0)+f( )2= , t= 2,h2=b a2= , (h2=h1/2)f(0) = ,f( ) = ,f( ) = ,1= (x)dx+ (x)dx= f(0)+f( )2+ f( )+f( )2= , t= 4,h3=b a4= , (h3=h2/2)17I3,1= [f(0) + 2f( ) + 2f( ) + 2f( ) +f( )] = , t= ,2=4I2,1 I1,14 1=43I2,1 13I1,1=43 13 = , t= ,2=4I3,1 I2,14 1=43I3,1 13I2,1=43 13 = , t= ,3=42I2,2 I1,242 1=1615I2,2 115I1,2=1615 115 = , t= 0 Figure 3: Example of Romberg integration18 PART II: Numerical DifferentiationFinite Divided Difference5 First Order Derivatives.

8 The first forward finite divided differenceUsing Taylor series,f(xi+1) =f(xi) +f (xi)h+f (xi)2!h2+O(h3)whereh=xi+1 xi. Thenf (xi)can be found asf (xi) =f(xi+1) f(xi)h+O(h)The first forward finite divided difference isf (xi) f(xi+1) f(xi)h The first backward finite divided differenceUsing Taylor series,f(xi 1) =f(xi) f (xi)h+f (xi)2!h2+O(h3)19whereh=xi xi 1. Thenf (xi)can be found asf (xi) =f(xi) f(xi 1)h+O(h)andf (xi)can also be approximated asf (xi) f(xi) f(xi 1)hwhich is called the first backward finite divided difference. The first centered finite divided differencef(xi+1) f(xi 1) = 2f (xi)h+O(h3)andf (xi)can be found asf (xi) =f(xi+1) f(xi 1)2h O(h2)andf (xi)can also be approximated asf (xi) f(xi+1) f(xi 1)2hwhich is called the first centered finite divided that the truncation error is of the order ofh2in contrast to the forwardand backward approximations that are of the order ofh.

9 Therefore, the centereddifference is a more accurate representation of the depiction of (a) forward, (b) backward, and (c) centered finite-divided-difference approximations of the first derivativeExample:Estimate the first derivative off(x) = + a step sizeh= Repeat the computation usingh= :The problem can be solved analyticallyf (x) = x ( ) = ,xi 1=xi h= 0, andf(xi 1) = ;xi= ,f(xi) = ;xi+1=xi+h= 1, andf(xi+1) = forward divided difference:f ( ) f(xi+1) f(xi)xi+1 xi= percentage relative error:| t|= ( ) ( ) 100% = backward divided difference:f ( ) f(xi) f(xi 1)xi xi 1= percentage relative error:| t|= ( ) ( ) 100% = centered divided difference:f ( ) f(xi+1) f(xi 1)xi+1 xi 1= percentage relative error:| t|= ( ) ( ) 100% = ,xi 1=xi h= , andf(xi 1) = ;xi= ,f(xi) =.

10 Xi+1=xi+h= , andf(xi+1) = forward divided difference:f ( ) f(xi+1) f(xi)xi+1 xi= percentage relative error:| t|= ( ) ( ) 100% = backward divided difference:f ( ) f(xi) f(xi 1)xi xi 1= percentage relative error:| t|= ( ) ( ) 100% = centered divided difference:f ( ) f(xi+1) f(xi 1)xi+1 xi 1= percentage relative error:| t|= ( ) ( ) 100% = centered finite divided difference and small step size achieves lower ap-proximation Order Derivatives: The second forward finite divided differencef(xi+2) =f(xi) +f (xi)(2h) +f (xi)2!(2h)2+O(h3)(2)f(xi+1) =f(xi) +f (xi)h+f (xi)2!h2+O(h3)(3)(2)-(3) 2:f(xi+2) 2f(xi+1) = f(xi) +f (xi)h2+O(h3)f (xi)can be found asf (xi) =f(xi+2) 2f(xi+1) +f(xi)h2+O(h)andf (xi)can be approximated asf (xi) f(xi+2) 2f(xi+1) +f(xi)h224 This is the second forward finite divided difference.