CH-04 (Pao - Engineering Analysis)

4Finite Differences,Interpolation, andNumerical Differentiation4.1 INTRODUCTIONLinear interpolation is discussed in the preceding chapter as a method for finding aparticular root of a polynomial, or, transcendental equation when the upper andlower bounds of the interval for search are provided. To continue the discussion ofthe general topic of interpolations which not necessarily linear could be quadratic(parabolic, cubic, quartic, and so on, we in this chapter present methods for thisgeneral need of interpolation in engineering analyses by treating not only equationsbut also a set of N tabulated data, (xi,yi) for i = 1–N. Finite difference table will beintroduced and constructed for the equally-spaced data, that is x2–x1 = x3–x2 = … =xN-xN–1.

This table can be utilized as a forward-difference, backward- difference, or,central-difference table depending on how its is applied for the interpolation use.Taylor’s series and a shifting operator are to be used in derivation of theinterpolation formulas in terms of the forward-difference, backward-difference, andcentral-difference operators. A program DiffTabl has been developed for printingout a difference table of a set of equally-spaced data.Differentiation operator will also be introduced for the derivation of the numerical differentiation needs. When a set of equally-spaced data, (xi,yi) for i = 1–N, aregiven, formulas in terms of the forward-difference, backward-difference, and centraldifference operators are derived for the need of calculating the value of dy/dx at alisted x value or unlisted. If x is not equal to one of the xi, interpolation anddifferentiation have to be done combinedly through a modification of the Taylor’sseries expansion.For curve-fit by polynomials and for interpolation, applications of the versatileLagrangian interpolation formula are also discussed.

A program called LagrangIis made available for this need.QuickBASIC, FORTRAN, and MATLAB versions of the above-mentionedprograms are to be provided. Application of the Mathematica’s function Interpolating Polynomial in place of LagrangI is demonstrated.In solution of the problems governed by a system of ordinary differential equations with either some initial and/or boundary conditions specified, the finite differences will be applied. In Chapter 6, such method for finding the approximate answerto the problem is discussed. Accuracy of such approximate solution will depend onthe increment of the independent variable, stepsize, adopted and on which approximate method is employed.© 2001 by CRC Press LLCBecause numerical differentiation is highly inaccurate, whenever possiblenumerical integration should be preferred over numerical differentiation.

In case thatone needs to find the velocity of a certain motion study and has the option ofcollecting the displacement or acceleration data, then the acceleration data shouldbe taken not the displacement data. The reason is that one has the choice of applyingnumerical differentiation to the displacement data or numerical integration to theacceleration data to obtain the velocity results. The numerical integration which isthe topic of Chapter 5 has the smoothing effect and hence is more accurate! Graphically, differentiation is of a local evaluation of determining the slope at a selectedpoint on a curve which could be the result of fitting a number of data points discussedin Chapter 3 while integration is of a global evaluation of finding the area under thecurve between two specified limits of the independent variable.

For a set of threegiven points fitted linearly by two linear segments and quadratically by a parabola,the slope at the mid-point could have very different slope values while the areasunder the linear segments and under the parabola would not differ too significantly.Hence, it is worthy of emphasizing that learning the computational methods is easierwhen compared to making decision of which method is best to solve the problemat hand.4.2 PROGRAM DIFFTABL — APPLICATIONSOF FINITE-DIFFERENCE TABLEProgram DiffTabl has been developed for the need of constructing a table of finitedifferences of a given set of N two-dimensional points, (xi,yi) for i = 1–N.

The xvalues are assumed to be equally spaced, i.e., , x2–x1 = x3–x2 = ••• = xN-xN–1 = h, hbeing called the increment, or, stepsize. This so-called difference table can be appliedfor interpolation of the y value for a specified, unlisted x value inside the range ofx = x1 and x = xN (extrapolation if outside the range), and differentiation. Table 1shows a typical difference table.The symbol used in Table 1 is called Forward Difference Operator. If we referthe numbers listed in the x and y columns as x1 to x6 and y1 to y6, respectively, thefirst number listed under y, 1.9495, is obtained from the calculation of y2–y1 andis identified as ∆y1. The last number listed in the y column, 5.3015, is equal toy6–y5 and referred to as y5.

Or, we may write the general formula as, for i = 1 to 5,∆y i = y i +1 − y i(1)yi is called the first forward difference of y at xi. The higher order forwarddifferences listed in Table 1 are obtained by extended application of Equation 1.That is,∆2 y i = ∆(y i +1 − y i ) = ∆y i +1 − ∆y i = y i +2 − 2 y i +1 + y i© 2001 by CRC Press LLC(2)TABLE 1Difference Table (y = 1 to 2x + 3x2 to 4x3 + 5x4).xy1.14.42651.26.3760y2y3y 4y 5y1.94950.6372.58651.38.96250.1260.7633.34951.412.31201.516.56251.621.86400.0120.1380.9014.25050.0000.0120.1501.0515.3015∆3y i = ∆2 (y i +1 − y i ) = ∆2 y i +1 − ∆2 y i= (y i +3 − 2 y i +2 + y i +1 ) − (y i +2 − 2 y i +1 + y i )(3)= y i +3 − 3y i +2 + 3y i +1 − y iand so on.

We shall show later how the third through seven columns of Table 1 canbe interpreted differently when the backward and central difference operators areintroduced. First, we will demonstrate how Table 1 can be applied for interpolationof the y value at an unlisted x value, say y(x = 1.24). To do that, the shifting operator,E, needs to be introduced. The definition of E is such that:Ey i = y i +1(4)That is, if E is operating on yi, the y value is shifted down to the next providedy value. Interpolation is a problem of not shifting a full step but a fractional step.For the need of finding y at x = 1.24, the x value falls between x2 = 1.2 and x3 =1.3. Since the stepsize, h, is equal to 0.1, a full shift from y2 = 6.3760 would leadto y3 which is equal to 8.9625.

We expect the value of y(x = 1.24) to be between y2and y3. Instead of E1y2, the value of E0.4y2 is to be calculated by shifting only 40%.To find the meaning of E0.24, or, more generally Er for 0<r<1, we substituteEquation 4 into Equation 1 to obtain:∆y i = y i +1 − y i = Ey i − y i = ( E − 1)y1∆ = E − 1,© 2001 by CRC Press LLCor,E = 1+ ∆(5,6)By application of binomial expansion, we can then have:E r = (1 + ∆ ) =r∞∑ ( )∆rkk(7)k =0where the binomial coefficients are defined as:( ) =1r0( )=rkandr( r − 1)…[ r − ( k − 1)]1⋅ 2L k(8)We can now use Equation 7 to obtain:y(x = 1.24) = E 0.4 y(x = 1.2) = E 0.4 y 2 = (1 + ∆ ) y 20.40.4(0.4 − 1) 2= 1 + 0.4 ∆ +∆ + … y 21⋅ 2((9))= 1 + 0.4 ∆ − 1.12 ∆2 + 0.064 ∆3 − 0.0416 ∆4 + 0.022952 ∆5 y 2Equation 9 can be applied for linear interpolation if up to the y2 terms areadopted; for parabolic interpolation if up to the 2y2 terms are adopted; and so on.Since Table 1 has up to the fifth order forward differences available but the lastcolumn contains a zero value, Equation 9 can therefore be effectively up to thefourth-order forward-difference interpolation.

The numerical results of y(x = 1.24)using linear, parabolic, cubic, and fourth-order are 7.4106, 7.3190, 7.3279, and7.3274, respectively. Since we know y = 1–2x + 3x2–4x3 + 5x4, the exact value ofy(x = 1.24) is equal to 7.3274.An explanation for discrepancies in all of these four attempts of interpolations,relative to the exact value, is provided in a homework exercise given in the Problemsset.BACKWARD-DIFFERENCE OPERATORNotice that the first numbers listed in columns three through seven in Table 1are the five forward differences of y1, and that only four forward differences (thesecond numbers in columns three through six) of y2 are available.

Lesser and lesserforward differences are available for later y’s until there is only y5 for y5. That isto say, to interpolate y(x) for an x value between x5 = 1.5 and x6 = 1.6, Equation 9can only be used up to the y5 term. To remedy this situation and to make most useof the provided set of 6 (x,y) data, it is appropriate at this time to introduce thebackward-difference operator, , which is defined as:∇y i = y i − y i −1© 2001 by CRC Press LLC(10)By combining Equations 1, 7, and 10, we notice that:∇y i +1 = y i +1 − y i = ∆y i(11)and()∇y i −1 = y i +1 − y i = 1 − E −1 y i +1(12)So,∇ = 1 − E,−1or,E −1 = 1 − ∇(13,14)Equation 12 is an important result because it indicates that the last numberslisted in columns three through seven of Table 1 are the first through five backwarddifferences of y6.

If we could derive an interpolation formula in terms of , thereare up to fifth-order backward difference of y6 available. Toward that end, let usconsider the need of interpolating the value of y(x = 1.56). This y value can bereached by shifting backward by 0.4 step from x = 1.6 since the stepsize for Table 1is h = 0.1. By using Equation 14 and noticing Equations 7 and 8, we can have:y(x = 1.56) = E −0.4 y(x = 1.6) = E −0.4 y 6 = (1 − ∇) y 60.40.4(0.4 − 1)= 1 + 0.4( −∇) +(−∇)2 + …y61⋅ 2()= 1 − 0.4∇ − 0.12∇2 − 0.064∇3 − 0.0416∇ 4 − 0.02295∇5 y 6 (15)One can then apply Equation 15 to obtain the interpolated y(x = 1.56) valuesusing up to the fifth order backward differences.