CH-02 (Pao - Engineering Analysis)

2Exact, Least-Squares, andCubic Spline Curve-Fits2.1 INTRODUCTIONEngineers conduct experiments and collect data in the laboratories. To make use ofthe collected data, these data often need to be fitted with some particularly selectedcurves. For example, one may want to find a parabolic equation y = c1 + c2x + c3x2which passes three given points (xi,yi) for i = 1,2,3. This is a problem of exact curvefit. Or, since knowing in advance that these three points should all fall on a straightline, but the reason that they are not is because of bad calibration of the measuringequipment or because of presence of noises in the testing environment.In case that we may want express this straight line by the equation y = c1 + c2xfor the stress and strain data collected for a stretching test of a metal bar in theelastic range, then the question of how to determine the two coefficients c1 and c2is a matter of deciding on which criterion to adopt.

The Least-Squares method isone of the criteria which is most popularly used. The two cases cited are theconsideration of adopting the two and three lowest polynomial terms, x0, x1, and x2,and linearly combining them.If the collected data are supposed to represent a sinusoidal function of time, thecurve to be determined may have to be assumed as x(t) = c1sint + c2sin3t + c3sin5t+ c4sin7t by linearly combining 4 odd sine terms.

This is the case of selecting fourparticular functions, namely, fi(t) = sin(2i–1)t for i = 1,2,3,4., and to determine thecoefficients c1–4 by application of the least-squares method.Often some special form of curve needs to be selected to fit a given set of data,the least-squares criterion can still be applied if mathematical transformations canbe found to convert the equation describing the curve into linear equations. This isdiscussed in a section devoted to transformed least-squares curve-fit.Another commonly applied curve-fit technique is the cubic spline method whichallows smooth cubic equations to be derived to ensure continuous slopes and curvatures passing all given points.

The mathematics involved in this method will bepresented.In the following sections, we shall discuss the development of the programsExactFit, LeastSq1, LeastSqG, and CubeSpln for the four curve-fit needs mentioned above.2.2 EXACT CURVE FITAs another example of solving a matrix equation, let us consider the problemof finding a parabolic equation y = c1 + c2x + c3x2 which passes three given points© 2001 by CRC Press LLC(xi,yi) for i = 1,2,3. This is a problem of exact curve-fit. By simple substitutions ofthe three points into the parabolic equation, we can obtain:c1 + c2 x i + c3x 2i = y ifor i = 1, 2, 3(1)In matrix form, we write these equations as:[A]{C} = {Y}(2)where {C} = [c1 c2 c3]T, {Y} = [y1 y2 y3]T, and [A] is a three-by-three coefficientmatrix whose elements if denoted as ai,j are to be calculated using the formula:a i, j = x ij−1for i, j = 1, 2, 3(3)It is easy to extend the above argument for the general case of exactly fitting Ngiven points by a N-1st degree polynomial y = c1 + c2x + ••• + cNxN–1.

The onlymodification needed is to allow the indices i and j in Equations 1 and 3 to be extendedto reach a value of N. A program called ExactFit has been prepared for this needby utilizing the subroutine GauJor to solve for the vector {C} from Equation 1 forthe general case of exactly fitting N given points. Listings for both FORTRAN andQuickBASIC versions along with sample numerical results are presented below.FORTRAN VERSION© 2001 by CRC Press LLCSample Applications© 2001 by CRC Press LLCQUICKBASIC VERSIONSample ApplicationMATLAB ApplicationFor handling the exact curve fit of N given points with a N-1st degree polynomial,there is no need to convert either the FORTRAN or QuickBASIC program ExactFit.

The sample problems therein can be readily solved by MATLAB as follows:© 2001 by CRC Press LLCNotice that the coefficient {C} for the curve-fit polynomial is obtained bysolving [A]{C} = {Y} where matrix [A] is formed by substituting the X valuesinto the xi terms for i = 0 to N–1 where N is the number of points provided.MATLAB function ones has been used to generate the first column of [A] andMATLAB matrix operation of C = A\Y which premultiplies {Y} by [A]–1 toobtain {C}.Also, this exact curve-fit problem can be treated as a special case of fitting Ngiven points by a linear combination of N selected functions which in this casehappens to be the polynomial terms of x0 to xN–1, by the least-squares method.

A mfile called LeastSqG.m which is discussed in the program LeastSqG can be readilyapplied to treat such a exact curve-fit problem. Here, we demonstrate how LeastSqG.m is used by MATLAB interactive operation in solving the sample problemspreviously presented in the FORTRAN and QuickBASIC versions of the programExactFit. First, a function must be prepared to describe the ith selected function xias follows:The results of four MATLAB applications are:© 2001 by CRC Press LLCNotice the coefficient vector {C} in the curve-fit polynomial p(x) = c1 + c2x +… + cNxN–1 is solved from the matrix equation [A]{C} = {R} where {A} and {R}are generated using the specified points based on the least squares criterion. Thesolution of [A]{C} = {R} is simply implemented by MATLAB with the statementC = A\R in the file LeastSqG.m.To verify whether the points have really been fitted exactly, Figure 1 is presented.It is plotted with the following MATLAB statements, adding to those already listedabove:Notice that for application of polyval.m, MATLAB needs the coefficients of thepolynomial arranged in descending order.

Since the array C contains the coefficients© 2001 by CRC Press LLCFIGURE 1. The parabolic curve passes through all of the three given points.in ascending order, a new array called Creverse is thus created for calculation of thecurve values for 1≤X≤3 and with an increment of 0.1. Figure 1 shows that theparabolic curve passes through all of the three given points.2.3 PROGRAM LEASTSQ1 — LEAST-SQUARES LINEAR CURVE-FITProgram LeastSq1 is designed for curve-fitting of a set of given data by a linearequation based on the least-squares criterion.2 If only two points are specified, alinear equation which geometrically describes a straight line can be uniquely derivedbecause the line must pass the two specified points.

This is the case of exact fit. (Seeprograms Gauss and LagrangI for examples of exact fit.) However, the specifieddata are often recorded from a certain experiment and due to inaccurate calibrationof equipment or due to environmental disturbances such as noise, heat, and so on,these data not necessarily follow an expected behavior which may be described bya type of predetermined equation. Under such a circumstance, a so-called forced fitis then required. As a simple example, supposing that we expect the measured setof three data points (Xi,Yi) for i = 1,2,3 to satisfy a linear law Y = c1 + c2X. If thesethree points happen to fall on a straight line, then we have a case of exact fit and© 2001 by CRC Press LLCthe coefficients c1 and c2 can be uniquely computed.

If these three points are not allon a straight line and they still need to be fitted by a linear equation Y = c1 + c2X,then they are forced to be fitted by a particular straight line based on a selectedcriterion, permitting errors to exist between the data points and the line.The least-squares curve-fit for finding a linear equation Y = c1 + c2X best representing the set of N given points, denoted as (Xi,Yi) for i = 1 to N, is to minimizethe errors between the straight line and the points to be as small as possible. Sincethe given points may be above, on, or below the straight line, and these errors maycancel from each other if they are added together, we consider the sum of the squaresof these errors. Let us denote yi to be the value of Y at Xi and S be the sum of theerrors denoted as ei for i = 1 to N, then we can write:NS = e12 + e 22 + … + e 2N =∑e2i(1)i =1where for i = 1,2,…,Ne i = y i − y i = c1 + c2 X i − Yi(2)It is obvious that since Xi and Yi are constants, the sum S of the errors is afunction of ci and c2.

To find a particular set of values of c1 and c2 such that S reachesa minimum, we may therefore base on calculus3 and utilize the conditions ∂S/∂c1 =0 and ∂S/∂c2 = 0. From Equation 1, we have the partial derivatives of S as: ∂e∂e ∂e∂S= 2 e1 1 + e 2 2 + … + e N N  = 2∂c1∂c∂c∂c1 11N∂e i∑ e ∂ci1i =1and ∂e∂e ∂e∂S= 2 e1 1 + e 2 2 + … + e N N  = 2∂c∂c∂c2∂c2 22N∂e i∑ e ∂cii =12From Equation 2, we note that ∂ei/∂c1 = 1 and ∂ei/∂c2 = Xi. Consequently, thetwo extremum conditions lead to two algebraic equations for solving the coefficientsc1 and c2:© 2001 by CRC Press LLCN∑i =11 c1 + N∑i =1X i  c2 =N∑Yii =1(3)andN∑i =1X i  c1 + N∑i =1X 2i  c2 =N∑X Yii(4)i =1Program LeastSq1 provided in both FORTRAN and QuickBASIC versions isdeveloped using the above two equations.

It can be readily applied for calculatingthe coefficients c1 and c2. Two versions are listed and sample applications arepresented below.QUICKBASIC VERSION© 2001 by CRC Press LLCSample ApplicationFORTRAN VERSION© 2001 by CRC Press LLCSample ApplicationMATLAB APPLICATIONA m file in MATLAB called polyfit.m can be applied to fit a set of given points(Xi,Yi) for i = 1 to N by a linear equation Y = C1X + C2 based on the least-squarescriterion. The function polyfit has three arguments, the first and second argumentsare the X and Y coordinate arrays of the given points, and the third argument specifiesto what degree the fitted polynomial is required. For linear fit, the third argumentshould be set equal to 1. The following shows how the results obtained for the sampleproblem used in the FORTRAN and QuickBASIC program LeastSq1:>> X = [1,2,3,5,8]; Y = [2,5,8,11,24]; A = polyfit(X,Y,1)C = 3.0195 – 1.4740If the third argument for the function polyfit is changed (from 1) to 2, 3, and4, we also can obtain the least-squares fits of the five given points with a quadratic,cubic, and quartic polynomials, respectively.