CH-04 (Pao - Engineering Analysis), страница 3

Input[3] stores the y array into the first column ofthe matrix yanddys by application of the Mathematica command Do. Such loopingis extended in Input[4] where the higher order differences are generated by usingan inner index i and an outer index j. The column number j of the matrix yanddysis increased from 2 to n but the length of each column is continuously decreased ton-j + 1. Such DoDo arrangement is made possible by keeping the y values and theirdifferences in a column-by-column form.4.3 PROGRAM LAGRANGI — APPLICATIONS OF LAGRANGIANINTERPOLATION FORMULAProgram LagrangI is designed to curve-fit a given set of n points, (xi,fi) for i = 1,2,…,n, by a polynomial of n-1st degree based on the Lagrangian InterpolationFormula:nf (x) =fi (x i − x k ) n∑ ∏[(x − x )i =1]kk =1k ≠i(1)If only the value of the function f(x) at a specified value of x = xs is needed,then Equation 1 can be applied to computef (xs ) =nfi n∑ ∏[(xi =1k =1k ≠i2(x i − x k ) ]− xk )(2)In Equations 1 and 2, the symbol is to represent a product of a specifiednumber of factors such as:n∏ F = F F …Fk1 2n(3)k =1Equation 1 can be proven if we write the equation which fits the n given points(xi,fi) for i = 1 to n by a combination of n functions L1 to n(x) as:f (x) = f1L1 (x) + f2 L 2 (x) + … + fn L n (x)(4)Notice that the ordinates f1 to n are utilized in Equation 4.

We expect the functionsL1 to n(x) to behave in such a way that when x = xi only the fiLi(x) term in Equation4 will contribute to f(x). That is to say when x = xi, Li(xi) should be equal to unityand the other L(x) should be equal to zero. Mathematically, we write demand that:© 2001 by CRC Press LLCL i (x i ) = 1L j (x i ) = 0 for j ≠ iand(5)The second condition of Equation 5 suggests that x-xk are factors of Li(x) fork = 1,2,…,n but not x-xi. Therefore, we may write:L i (x) = ci (x − x1 )(x − x 2 )…(x − x i −1 )(x − x i +1 )…(x − x n )(6)The constant associated with Li(x), ci is to be determined by satisfying the firstcondition of Equation 5.

That is:ci = 1 n∏ (x − x )i(7)kk =1k ≠iConsequently, the complete expression for Li(x) is:∏[(x − x ) (x )]nL i (x) =k(8)i−xkk =1k ≠iAnd, when Equation 8 is substituted into Equation 4, we arrive at Equation 1.A numerical example will clarify the application of Equation 2. Consider thecase of three given points (x1,f1) = (1,2), (x2,f2) = (1.5,2.5), (x3,f3) = (3,4), then n =3. If we need to calculate f(x = 2), Equation 2 can be used to find the equation whichpasses all three points. That is:f (x) =(x − x2 )(x − x3 ) f + (x − x1 )(x − x3 ) f(x1 − x2 )(x1 − x3 ) 1 (x2 − x1 )(x2 − x3 ) 2+(x − x1 )(x − x2 ) fx( 3 − x1 )(x3 − x2 ) 3=(x − 1.5)(x − 3) 2 + (x − 1)(x − 3) 2.5(1 − 1.5)(1 − 3)(1.5 − 1)(1.5 − 3)+(x − 1)(x − 1.5) 4(3 − 1)(3 − 1.5)= 2(x − 4.5x + 4.5) −= x +1© 2001 by CRC Press LLC()()10 24x − 4 x + 3 + x 2 − 2.5x + 1.533When x = 2, f(x = 2) = 3.

Actually, the value of f(x = 2) can be specificallycalculated as:f (x = 2) =(x − x2 )(x − x3 ) f + (x − x1 )(x − x3 ) f(x1 − x2 )(x1 − x3 ) 1 (x2 − x1 )(x2 − x3 ) 2+(x − x1 )(x − x2 ) f(x3 − x1 )(x3 − x2 ) 3=(2 − 1.5)(2 − 3) 2 + (2 − 1)(2 − 3) 2.5(1 − 1.5)(1 − 3)(1.5 − 1)(1.5 − 3)+(2 − 1)(2 − 1.5) 4(3 − 1)(3 − 1.5)= −1 +QUICKBASIC VERSION© 2001 by CRC Press LLC10 2+ =33 3Sample ApplicationFORTRAN VERSION© 2001 by CRC Press LLCSample ApplicationMATLAB APPLICATIONA m file called LagrangI.m can be created and added to MATLAB m files forinterpolating a Y value for a giving X value based on a set of (X,Y) data point usingthe Lagrangian formula. This file may be written as:This m file can then be applied by specifying the data points, X vs.

Y, asillustrated by the following examples:The graphic capability of MATLABcan also be utilized here to interpretLagrangian interpolation. In Figure 6, five given points marked with the character* have been exactly fitted with a fourth-order polynomial which is plotted for 1≤X≤5with a solid line. The interpolation at X = 4.56 using Lagrangian formula is illustratedby the broken line and dotted line. The interactively entered MATLAB statements,in addition to those already displayed above, are:© 2001 by CRC Press LLCFIGURE 6. Five given points marked with the character * have been exactly fitted with afourth-order polynomial which is plotted for 1≤X≤5 with a solid line.

The interpolation atX = 4.56 using Lagrangian formula is illustrated by the broken line and dotted line.© 2001 by CRC Press LLCNotice that plot.m automatically uses solid, broken, and dotted lines to plot thefour-order polynomial curve based on arrays XC and YC, and the vertical line basedon arrays XV, YV, and the horizontal line based on arrays XH, YH, respectively.The details involved in exact curve-fit of the five given point by applying LeastSqG.m already has been discussed in the program Gauss.

The coefficients, {C}, ofthe fourth-order polynomial determined by LeastSqG.m are arranged in descendingorder. In order to apply polyval.m of MATLAB, the order of {C} has to be reversedand stored in {Creverse} which is implemented above by the for and end loop.MATHEMATICA APPLICATIONSDerivation of the polynomial which passing through a set of given (x,y) pointsbased on the Lagrangian formula can be achieved by application of the InterpolatingPolynomial function of Mathematica. For example, a fourth-order polynomial canbe derived for a given set of 5 (x,y) data points as follows:In[1]: = pofx = InterpolatingPolynomial [{{1,2},{2,4},{3,6},{4,8},{5,11}},x](−4 + x) (−3 + x) (−2 + x)  −1 + xOut[1]: = 2 +  2 +)(24To interpolate the y value of using the derived polynomial at x equal to 4.56,we replace all x’s appearing in the above expression (saved in pofx) with a valueof 4.56 by interactively enteringIn[2]: = pofx/.

x -> 4.56Out[2]: = 9.45174Linear and parabolic interpolations can also be implemented by selecting appropriate data points from the given set. For example, to interpolate the y value at x =1.25 by linear interpolation, we enter:In[3]: = p1 = InterpolatingPolynomial[{{1,2},{2,4}},x]Out[3]: = 2 + 2 (–1 + x)In[4]: = p1/. x -> 1.25Out[4]: = 2.5To parabolically interpolate the y value at x = 3.75 using the points (3,6), (4,8),and (5,11), the interactive application of Mathematica goes as:In[5]: = p2 = InterpolatingPolynomial[{{3,6},{4,8},{5,11}},x]© 2001 by CRC Press LLCOut[5]: = 6 +  2 +−4 + x (−3 + x)2 In[6]: = p2/. x -> 3.75 Out[6]: = 7.406254.4 PROBLEMSDIFFTABL1.

Construct the difference table based on the following listed data and thenfind the y value at x = 4.5 by using the backward-difference formula upto the third-order difference.xy1224374125202. Explain why interpolations using Equation 9 by the first through fourthorders all fail to match the exact value of y(x = 1.24) = 7.3274 by making4 plots for x values ranging from 1.2 to 1.3 with an increment of x =0.001. These 4 plots are to be generated with the 4 equations obtainedwhen the first 2, 3, 4, and 5 points are fitted by a first-, second-, third-,and fourth-degree polynomials, respectively. Also, draw a x = 1.24, vertical line crossing all 4 curves.3.

Find the first-, second-, third-, and fourth-order results of y(x = 1.56) byuse of Equation 15.4. Write Er in terms of binomial coefficient and the backward-differenceoperator , similar to Equation 7.5. Find the first-, second-, third-, and fourth-order results of y(x = 1.24) byuse of Equation 24.6.

Find the first-, second-, third-, and fourth-order results of y(x = 1.56) byuse of Equation 25.7. Given 6 (x,y) points (1,0.2), (2,0.4), (3,0.7), (4,1.5), (5,2.9), and (6,4.7),parabolically interpolate y(x = 3.4) first by use of forward differences andthen by use of backward differences.8. Modify either the QuickBASIC or FORTRAN version of the programDiffTabl to include the fifth difference for the need of forward or backward interpolation and numerical differentiation.9. Given 5 (x,y) points (0,0), (1,1), (2,8), (3,27), and (4,64), construct acomplete difference table based on these data. Compute (1) y value at x =1.25 using a forward, parabolic (second-order) interpolation, (2) y valueat x = 3.7 using a backward, cubic (third-order) interpolation, and (3)dy/dx value at x = 0 using a forward, third-order approximation.10.

Based on Equation 21, derive the forward-difference formulas for D2yiand D3yi.11. Use the result of Problem 10 to compute D2y2 and D3y1 by adopting theforward-difference terms in Table 1 as high as available.© 2001 by CRC Press LLC12. Use the data in Table 1 to compute the first derivative of y at x = 1.155by including terms up to the third-order forward difference.13.

Apply MATLAB for the points given in Problem 1 to print out the rowsof x, y, ∆y, 2y, 3y, and 4y.14. Same as Problem 13 but the points in Problem 6.15. Apply Mathematica and DO loops to print out a difference table similarto that shown in Mathematica Application of Section 4.2 for the pointsgiven in Problem 1.16.

Apply Mathematica and DO loops to print out a difference table similarto that shown in Mathematica Application of Section 4.2 for the pointsgiven in Problem 6.17. Compute the binomial coefficient for r = 0.4 and k = 1,2,3,4,5 accordingto Equation (8) in Section 4.2 using MATLAB.18. Rework Problem 17 but using Mathematica.LAGRANGI1. Given five points (1,1), (2,3), (3,2), (4,5), and (5,4), use the last threepoints and Lagrangian interpolation formula to compute y value at x = 6.2.