Просмотреть файл в отдельном окне: 097db74dac7bf547a1327c4e1338a012.pdf

2

Exact, Least-Squares, and

Cubic Spline Curve-Fits

2.1 INTRODUCTION

Engineers conduct experiments and collect data in the laboratories. To make use of

the collected data, these data often need to be fitted with some particularly selected

curves. For example, one may want to find a parabolic equation y = c1 + c2x + c3x2

which 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 straight

line, but the reason that they are not is because of bad calibration of the measuring

equipment 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 + c2x

for the stress and strain data collected for a stretching test of a metal bar in the

elastic range, then the question of how to determine the two coefficients c1 and c2

is a matter of deciding on which criterion to adopt. The Least-Squares method is

one of the criteria which is most popularly used. The two cases cited are the

consideration 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, the

curve 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 four

particular functions, namely, fi(t) = sin(2i–1)t for i = 1,2,3,4., and to determine the

coefficients 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 can

be found to convert the equation describing the curve into linear equations. This is

discussed in a section devoted to transformed least-squares curve-fit.

Another commonly applied curve-fit technique is the cubic spline method which

allows smooth cubic equations to be derived to ensure continuous slopes and curvatures passing all given points. The mathematics involved in this method will be

presented.

In the following sections, we shall discuss the development of the programs

ExactFit, LeastSq1, LeastSqG, and CubeSpln for the four curve-fit needs mentioned above.

2.2 EXACT CURVE FIT

As another example of solving a matrix equation, let us consider the problem

of 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 of

the three points into the parabolic equation, we can obtain:

c1 + c2 x i + c3x 2i = y i

for 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 coefficient

matrix whose elements if denoted as ai,j are to be calculated using the formula:

a i, j = x ij?1

for i, j = 1, 2, 3

(3)

It is easy to extend the above argument for the general case of exactly fitting N

given points by a N-1st degree polynomial y = c1 + c2x + ••• + cNxN–1. The only

modification needed is to allow the indices i and j in Equations 1 and 3 to be extended

to reach a value of N. A program called ExactFit has been prepared for this need

by utilizing the subroutine GauJor to solve for the vector {C} from Equation 1 for

the general case of exactly fitting N given points. Listings for both FORTRAN and

QuickBASIC versions along with sample numerical results are presented below.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

© 2001 by CRC Press LLC

QUICKBASIC VERSION

Sample Application

MATLAB Application

For 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 LLC

Notice that the coefficient {C} for the curve-fit polynomial is obtained by

solving [A]{C} = {Y} where matrix [A] is formed by substituting the X values

into 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] and

MATLAB matrix operation of C = A\Y which premultiplies {Y} by [A]–1 to

obtain {C}.

Also, this exact curve-fit problem can be treated as a special case of fitting N

given points by a linear combination of N selected functions which in this case

happens to be the polynomial terms of x0 to xN–1, by the least-squares method. A m

file called LeastSqG.m which is discussed in the program LeastSqG can be readily

applied to treat such a exact curve-fit problem. Here, we demonstrate how LeastSqG.m is used by MATLAB interactive operation in solving the sample problems

previously presented in the FORTRAN and QuickBASIC versions of the program

ExactFit. First, a function must be prepared to describe the ith selected function xi

as follows:

The results of four MATLAB applications are:

© 2001 by CRC Press LLC

Notice 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. The

solution of [A]{C} = {R} is simply implemented by MATLAB with the statement

C = 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 listed

above:

Notice that for application of polyval.m, MATLAB needs the coefficients of the

polynomial arranged in descending order. Since the array C contains the coefficients

© 2001 by CRC Press LLC

FIGURE 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 the

curve values for 1?X?3 and with an increment of 0.1. Figure 1 shows that the

parabolic curve passes through all of the three given points.

2.3 PROGRAM LEASTSQ1 — LEAST-SQUARES LINEAR CURVE-FIT

Program LeastSq1 is designed for curve-fitting of a set of given data by a linear

equation based on the least-squares criterion.2 If only two points are specified, a

linear equation which geometrically describes a straight line can be uniquely derived

because the line must pass the two specified points. This is the case of exact fit. (See

programs Gauss and LagrangI for examples of exact fit.) However, the specified

data are often recorded from a certain experiment and due to inaccurate calibration

of 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 by

a type of predetermined equation. Under such a circumstance, a so-called forced fit

is then required. As a simple example, supposing that we expect the measured set

of three data points (Xi,Yi) for i = 1,2,3 to satisfy a linear law Y = c1 + c2X. If these

three points happen to fall on a straight line, then we have a case of exact fit and

© 2001 by CRC Press LLC

the coefficients c1 and c2 can be uniquely computed. If these three points are not all

on 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 selected

criterion, 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 minimize

the errors between the straight line and the points to be as small as possible. Since

the given points may be above, on, or below the straight line, and these errors may

cancel from each other if they are added together, we consider the sum of the squares

of these errors. Let us denote yi to be the value of Y at Xi and S be the sum of the

errors denoted as ei for i = 1 to N, then we can write:

N

S = e12 + e 22 + … + e 2N =

?e

2

i

(1)

i =1

where for i = 1,2,…,N

e 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 a

function of ci and c2. To find a particular set of values of c1 and c2 such that S reaches

a 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 ?

?

1

1

N

?e i

? e ?c

i

1

i =1

and

? ?e

?e ?

?e

?S

= 2? e1 1 + e 2 2 + … + e N N ? = 2

?

c

?

c

?c2

?c2 ?

?

2

2

N

?e i

? e ?c

i

i =1

2

From Equation 2, we note that ?ei/?c1 = 1 and ?ei/?c2 = Xi. Consequently, the

two extremum conditions lead to two algebraic equations for solving the coefficients

c1 and c2:

?

?

?

© 2001 by CRC Press LLC

N

?

i =1

?

?

1? c1 + ?

?

?

N

?

i =1

?

X i ? c2 =

?

N

?Y

i

i =1

(3)

and

?

?

?

N

?

i =1

?

?

X i ? c1 + ?

?

?

N

?

i =1

?

X 2i ? c2 =

?

N

?X Y

i

i

(4)

i =1

Program LeastSq1 provided in both FORTRAN and QuickBASIC versions is

developed using the above two equations. It can be readily applied for calculating

the coefficients c1 and c2. Two versions are listed and sample applications are

presented below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

MATLAB APPLICATION

A 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-squares

criterion. The function polyfit has three arguments, the first and second arguments

are the X and Y coordinate arrays of the given points, and the third argument specifies

to what degree the fitted polynomial is required. For linear fit, the third argument

should be set equal to 1. The following shows how the results obtained for the sample

problem 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.4740

If the third argument for the function polyfit is changed (from 1) to 2, 3, and

4, we also can obtain the least-squares fits of the five given points with a quadratic,

cubic, and quartic polynomials, respectively. When the third argument is set equal

to 4, we then have the case of exact curve-fit of five points by a fourth-order

polynomial. Readers are referred to the program ExactFit for more discussions.

Also, it is of interest to know whether one may select an arbitrary set of

functions and linearly combine them for least-squares fit, instead of the unbroken

set of polynomial terms X0, X1, X2, …, XN. Program LeastSqG to be presented

in the next section will discuss such generalized least-squares curve-fit. But before

we do that, let us first look into a situation where program LeastLQ1 can be

applied for a given set of data after some mathematical transformations are

employed to modify the data.

Transformed Least-Squares Curve-Fit

There are occasions when we know in advance that a given set of data supposed

to fall on a curve described by exponential equations of the type:

© 2001 by CRC Press LLC

y = b1e b2x

(5)

y = c1xe c2x

(6)

or

To determine the coefficients b1 and b2, or, c1 and c2 based on the least-squares

criterion, Equation 5 or 6 need to be first transformed into a linear form. To do so,

let us first consider Equation 5 and take natural logarithm of both sides. It gives:

n y = n b1 + n e b2x = n b1 + b 2 x

(7)

If we introduce new variable z, and new coefficients a1 and a2 such that:

z = n y,

a1 = n b1,

and

a 2 = b2

(8,9,10)

Then Equation 7 becomes:

z = a1 + a 2 x

(11)

Hence, if we need to determine the coefficients b1 and b2 for Equation 2.8 based

on N pairs of (xi,yi), for i = 1 to N, values and on the least-squares criterion, we

simply generate N z values according to Equation 2.11 and then use the N (xi,zi)

values as input for the program LeastSq1 and expect the program to calculate a1 and

a2. Equations 9 and 10 suggest that b2 is to have the value of a2 while b1 should be

equal to e raised to the a1 power, or, EXP(a1) with EXP being the exponential function

available in the FORTRAN or QuickBASIC library.

Equation 6 can be treated in a similar manner by taking logarithms of both sides

to obtain:

n y = n c1 + n x + n e c2x = n c1 + n x + c2 x

or

n y? n x = n

y

= n c1 + c2 x

x

(12)

If we introduce new variable z, and new coefficients a1 and a2 such that:

z = n(y x),

a1 = n c1,

and

a 2 = c2

(13,14,15)

Then, Equation 12 becomes Equation 11 and a1 and a2 can be obtained by the

program LeastSq1 using the data set of (xi,yi/xi) values.

As a numerical example, consider the case of a set of nine stress-versus-strain

( vs. ) data collected from a stretching test of a long bar: (.265,1025), (.4,1400),

(.5,1710), (.7,2080), (.95,2425), (1.36,2760), (2.08,3005), (2.45,2850), and

© 2001 by CRC Press LLC

(2.94,2675) where the units for the strains and stresses are in microinch x 102 and

lb/in2, respectively. When program LeastSq1 is applied for the modified data of

(,), the coefficients for the linear fit are a1 = 15.288 and a2 = –537.71. Consequently, according to Equations 13 and 14, and realizing that x and y are now and

, respectively, we have arrived at = 4.3615 ? 106e–537.71. The derived curve and

the given points are plotted in Figure 2 which shows the curve passes the origin as

it should.

FIGURE 2. The curve passes the origin as it should.

© 2001 by CRC Press LLC

2.4 PROGRAM LEASTSQG — GENERALIZED

LEAST-SQUARES CURVEFIT

Program LeastSqG is designed for curve-fitting of a set of given data by a linear

combination of a set of selected functions based on the least-squares criterion.2

Let us consider N points whose coordinates are (Xk,Yk) for k = 1 to N and let

the M selected function be f1(X) to fM(X) and the equation determined by the leastsquares curve-fit be:

M

Y( X) = a1f1 ( X) + a 2 f2 ( X) + … + a M fM ( X) =

? a f (X)

j j

(1)

j=1

The least-squares curve-fit for finding the coefficients c1–M is to minimize the

errors between the computed Y values based on Equation 1 and the given Y values

at all Xk’s for k = 1 to N. Let us denote yk to be the value of Y calculated at Xk

using Equation 1 and S be the sum of the errors denoted as ek for k = 1 to N. Since

the yk could either be greater or less than Yk, these errors ek’s may cancel from each

other if they are added together. We therefore consider the sum of the squares of

these errors and write:

N

S = e12 + e 22 + … + e 2N =

?e

2

k

(2)

k =1

where for k = 1,2,…,N

?

e k = y k ? Yk = ?

?

?

M

?

j=1

?

a jfj ( X k )? ? Yk

?

?

(3)

It is obvious that since Xk and Yk are constants, the sum S of the errors is a

function of a1 to M. To find a particular set of values of a1 to M such that S reaches a

minimum, we may therefore base on calculus3 and utilize the conditions ?S/?ai = 0

for i = 1 to M. From Equation 2, we can have the partial derivatives of S with respect

to ai’s, for i = 1 to M, as:

N

? ?e

?e ?

?e

?e

?S

= 2? e1 1 + e 2 2 + … + e N N ? = 2

ek k

?a i

?a i

?a i

?a i ?

? ?a1

k =1

?

From Equation 3, we note that ?ek/?ai = fi(Xk). Consequently, the M extremum

conditions, ?S/?ai = 0 for i = 1 to M, lead to M algebraic equations for solving the

coefficients a1 to M. That is, for i = 1 to M:

M

?

? N

fi ( X k )fj ( X k )?a j =

?

??

?? k =1

??

j=1

© 2001 by CRC Press LLC

N

? f (X )Y

i

k =1

k

k

(4)

If we express Equation 4 in matrix notation, it has the simple form:

[C]{A} = {R}

(5)

where [C] is a MxM square coefficient matrix, and {A} and {R} are Mx1 column

matrices (vectors). {A} contains the unknown coefficients a1 to M needed in Equation 1.

If we denote the elements in [C] and {R} as cij and ri, respectively, Equation 5

indicates that these elements are to be calculated using the formulas, for i = 1,2,…,M:

N

cij =

? f (X )f (X )

i

k

j

k

(6)

k =1

and

N

ri =

? f (X )Y

i

k

k

(7)

k =1

The above derivation appears to be too mathematical; a few examples of actual

curve-fit will clarify the procedure involved. As a first example, consider the case

of selecting two (M = 2) functions f1(X) = 1 and f2(X) = X to fit three given points

(N = 3), (X1,Y1) = (1,1), (X2,Y2) = (2.6,2), and (X3,Y3) = (2.8,2). Equations 6 and 7

then provide the following:

and

Hence, the system of two linear algebraic equations for finding a1 and a2 for the

straight-line equation is:

?3

?6.4

?

© 2001 by CRC Press LLC

6.4 ? ? a1 ? ? 5 ?

=

15.6?? ??a 2 ?? ??11.8??

The solution can be obtained by application of Cramer’s Rule to be a1 = 0.42466

and a2 = 0.58219. More examples will be given after we discuss how computer

programs can be written to compute [C] and {R} and then solve for {A}.

Program LeastSqG provided in both FORTRAN and QuickBASIC versions

is developed using the above two equations. It can be readily applied for calculating

the coefficients c1 to M. Both QucikBASIC and FORTRAN versions are listed and

sample applications are presented below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Applications

When four functions are selected as those listed in SUB FS, an interactive

application of the program LeastSqG QuickBASIC version using the input data

entered through keyboard has resulted in a screen display of:

If three sinusoidal functions, sin(x/20), sin(3x/20), and sin(5x/20) were

selected and replacing those listed in SUB FS, another interactive application of the

program LeastSqG QuickBASIC version is shown below.

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

By selecting the four functions listed in Subroutine FS, an interactive application

of the program LeastSqG using the input data given below has resulted in a screen

display of:

MATLAB APPLICATION

A LeastSqG.m file can be created and added to MATLAB m files which will

take N sets of X and Y points and fitted by a linear combination of M selected

functions in the least-squares sense. The selected functions can be specified in

another m file called FS.m (using the same name as in the FORTRAN and QuickBASIC versions). First, let us look at a version of LeastSqG.m:

Notice that the coefficients C’s is obtained by solving [A]{C} = {R} as in the

text. For MATLAB, a simple matrix multiplication of the inverse of [A] to and on

the left of the vector {R}. Complete execution of LeastSqG.m will be indicated by

a display of the matrix [A], vector {R}, and the solution vector {C}. The expression

feval(funfcn,i,X(k)) in the above program is to evaluate the ith function at X(k)

defined in a function file to be specified when LeastSqG.m is applied which is to

be illustrated next.

© 2001 by CRC Press LLC

Consider the case of given 5 (X,Y) points (N = 5) which are (1.4,2.25), (3.2,15),

(4.8,26.25), (8,33), and (10,35). And, the selected functions are sin(X/20),

sin(3X/20), and sin(5X/20). That is, M = 3. The supporting function FS.m is

simply:

function F=FS(i,xv)

F=sin((2.*i-1).*xv.(pi./20);

Having prepared this file FS.m on a disk which is in drive A, we can now apply

LeastSqG.m interactively as follows:

>> X = [1.4,3.2,4.8,8,10]; Y = [2.25,15,26.25,33,35]; C = LeastSqG(A:FS,X,Y,3,5)

The results shown on screen are:

If four functions X, X2, sin(X), and cos(X) are selected, we may change the

FS.m file to:

© 2001 by CRC Press LLC

Same as for the FORTRAN and QuickBASIC versions, if we are given six

(X,Y) points, (1,2), (2,4), (3,7), (4,11), (5,23), and (6,45), MATLAB application of

LeastSqG.m will be:

>> X=[1,2,3,4,5,6]; Y=[2,4,7,11,23,45]; C-LeastSqG('A:FS',X,Y,4,6)

The results are:

Notice that the values in [A] and {R} are to be multiplied by the factor 1.0e +

003 as indicated. For saving space, [A], {R}, and {C} are listed side-by-side when

actually they are displayed from top-to-bottom on screen. To further utilize the

capability of MATLAB, the obtained {C} values are checked to plot the fitted curve

against the provided six (X,Y) points. The following additional statements are

entered in order to have a screen display as illustrated in Figure 3:

FIGURE 3. The ‘*’ argument in the second plot statement requests that the character * is

to be used for plotting the given points,

© 2001 by CRC Press LLC

The statement XC = [1:0.2:6] generates a vector of XC containing value from

1 to 6 with an increment of 0.2. The command “hold” enables the first plot of XC

vs. YC which is the resulting least-squares fitted curve, to be held on screen until

the given points (X,Y) are superimposed. The ‘*’ argument in the second plot

statement requests that the character * is to be used for plotting the given points, as

illustrated in Figure 3.

Next, another example is given to show how MATLAB statements can be applied

directly with defining a m function, such as FS which describes the selected set of

functions for least-squares curve fit. Consider the problem of least-squares fit of the

points (1,2), (3,5), and (4,13) by the linear combination Y = C1f1(X) + C2f2(X) where

f1(X) = X–1 and f2(X) = X2. An interactive application of MATLAB may go as follows:

The resulting curve is plotted in Figure 4 using 31 points, (XC,YC), calculated

based on the equation Y = –7.0526(X–1) + 2.1316X2 for X values ranging from 1

to 4. The three given points, (X,Y), are superimposed on the graph using ‘*’

character. Notice from the above statements, the coefficients C(1) and C(2) are solved

from the matrix equation [A]{C} = {R} where the elements in [A] are generated

using interactively entered statement and so are the elements of {R}. MATLAB

matrix operations such as transposition, subtraction, multiplication, and inversion

are all involved. Also, notice that here no use of MATLAB ‘hold’ as for generating

Figure 1, is necessary if X, Y, XC, and YC are all used as arguments in calling the

© 2001 by CRC Press LLC

FIGURE 4. The curve is plotted using 31 points, (XC,YC), calculated based on the equation

Y = –7.0526(X–1) + 2.1316X2 for X values ranging from 1 to 4.

plot function. The statement XC = 1:0.1:4 generates a XC row matrix having elements starting from a value equal to 1, equally incremented by 0.1 until a value of

4 is reached.

MATHEMATICA APPLICATIONS

Mathematica has a function called Fit which least-squares fits a given set of

(x,y) points by a linear combination of a number of selected functions. As the first

example of interactive application, let us find the best straight line which gives the

least squared errors for fitting a set of five points. This is the case of two selected

functions f1(x) = 1 and f2(x) = x. the interactive application goes as follows:

In[1]: = Fit[{{1,2},{2,5},{3,8},{5,11},{8,24}}, {1, x}, x]

Out[1]: = –1.47403 + 3.01948 x

Notice that Fit has three arguments: first argument specifies the data set, second

argument lists the selected function, and the third argument is the variable for the

© 2001 by CRC Press LLC

derived equation. In case that three points are given and the two selected functions

are f1(x) = x–1 and f2(x) = x2, then the interactive operation goes as follows:

In[2]: = Fit[{{1,2},{3,5},{4,13}}, {x–1, x^2}, x]

Out[2]: = –7.05263 (–1 + x) + 2.13158 x

Two other examples previously presented in the MATLAB applications can also

be considered and the results are:

In[3]: = (Fit[{{1,2}, {2,4}, {3,7}, {4,11}, {5,23}, {6,45}},

{x, x^2, Sin[x], Cos[x]}, x])

Out[3]: = –4.75756 x + 2.11159 x2 – 0.986915 Cos[x] + 5.76573 Sin[x]

and

In[4]: = (Fit[{{1.4, 2.25}, {3.2, 15}, {4.8, 26.25}, {8, 33}, {10, 35}},

{Sin[Pi*x/20], Sin[3*Pi*x/20], Sin[5*Pi*x/20]}, x])

Pi x ?

3 Pi x ?

5 Pi x ?

Out[4]: = 35.9251 Sin ??

? 1.19261 Sin ??

? 3.46705 Sin ??

? 20 ??

? 20 ??

? 20 ??

All of the results obtained here are in agreement with those presented earlier.

2.5 PROGRAM CUBESPLN — CURVE FITTING

WITH CUBIC SPLINE

If a set of N given points (Xi,Yi) for i = 1, 2,…,N is to be fitted with a polynomial

of N–1 degree passing these points exactly, the polynomial curve will have many

fluctuations between data points as the number of points, N, increases. A popular

method for avoiding such over-fluctuation is to fit every two adjacent points with a

cubic equation and to impose continuity conditions at the point shared by two

neighboring cubic equations. This is like using a draftsman’s spline attempting to

pass through all of the N given points. It is known as cubic-spline curve fit. For a

typical interval of X, say Xi to Xi + 1, the cubic equation can be written as:

Y = a i X3 + bi X 2 + ci X + d i

(1)

where ai, bi, ci, and di are unknown coefficients. If there are N–1 intervals and each

interval is fitted by a cubic equation, the total number of unknown coefficients is

equal to 4(N–1). They are to be determined by the following conditions:

1. Each curve passes two specified points. The first and last specified points

are used once by the first and N-1st curves, respectively, whereas all

interior points are used twice by the two cubic curves meeting there. This

gives 2 + 2(N–2), or, 2N–2 conditions.

© 2001 by CRC Press LLC

2. Every two adjacent cubic curves should have equal slope (first derivative

with respect to X) at the interior point which they share. This gives N–2

conditions.

(2)

3. For further smoothness of the curve fit, every two adjacent cubic curves

should have equal curvature (second derivative with respect to X) at the

interior point which they share. This gives N–2 conditions.

4. At the end points, X1 and XN, the nature spline requires that the curvatures

be equal to zero. This gives 2 conditions.

Instead of solving the coefficients in Equation 1, the usual approach is to apply

the above-listed conditions for finding the curvatures, Y? at the interior points, that

is for i = 2,3,…,N–1 since (Y?)1 = (Y?)N = 0. To do so, we notice that if Y is a cubic

polynomial of X, Y? is then linear in X and can be expressed as:

Y?? = AX + B

(3)

If this is used to fit the ith interval, for which the increment in X is here denoted

as Hi = Xi + 1–Xi, we may replace the unknown coefficients A and B with the

curvatures at Xi and Xi + 1, (Y?)i and (Y?)i + 1 by solving the two equations:

Yi?? = AX i + B

and

Yi??+1 = AX i +1 + B

(4,5)

By using the Cramer’s rule, it is easy to obtain:

A=

Yi?+1 ? Yi??

X i +1 ? X i

and

B = Y?? ?

X i (Yi??+1 ? Yi??)

X i +1 ? X i

(6,7)

Consequently, Equation 3 can be written as:

Y?? =

Yi??

(X ? X) + YHi??+1 (X ? Xi )

H i i +1

i

(8)

Equation 8 can be integrated successively to obtain the expressions for Y and

Y as:

2

2

? Yi??

Y??

Y? +

X i +1 ? X) + i +1 ( X ? X i ) + C1

(9)

(

2Hi

2H i

and

3

3

Y??

Y??

Y = i ( X i +1 ? X) + i +1 ( X ? X i ) + C1X + C2

(10)

6H i

6H i

© 2001 by CRC Press LLC

The integration constants C1 and C2 can be determined by the conditions that at

Xi, Y = Yi and at Xi + 1, Y = Yi + 1. Based on Equation 10, the two conditions lead to:

Yi =

Yi?? 2

H + C1X i + C2

6 i

and

Yi +1 =

Yi??+1 2

H + C1X i +1 + C2

6 i

(11,12)

Again, Cramer’s rule can be applied to obtain:

C1 =

(Y??? Y?? )H

i +1

i

6

i

?

Yi ? Yi +1

Hi

(13)

and

C2 = ?

(Xi+1Yi??? XiYi??+1 )Hi ? XiYi+1 ? Xi+1Yi

6

Hi

Hi

(14)

Substituting C1 and C2 into Equations 11 and 12 and rearranging into terms

involving the unknown curvatures, the expressions for the cubic curve are:

Y=

3

? Y?? ? ( X ? X )3

?

Yi?? ? ( X i +1 ? X)

i

?

? H i ( X i +1 ? X)? + i +1 ?

? H i ( X ? X i )?

Hi

6 ?

6 ? Hi

??

??

?

?

? X ? X?

? X ? Xi ?

+ Yi ? i +1

+ Yi +1 ?

?

?

? Hi ?

? Hi ?

(15)

? H ( X ? X )2 ?

? ( X ? X )2 H ? Y ? Y

i

i

? + Yi??+1 ?

Y? = Yi?? ? i ? i +1

? i ? + i +1

Hi

2Hi

6?

??

?? 6

?? 2 H i

?

(16)

and

Equation 15 indicates that the cubic curve for the ith interval is completely

defined if in addition to the specified values of Yi and Yi + 1, the curvatures at Xi and

Xi + 1, (Y?)i and (Y?)i + 1 respectively, also can be found. To find all of the curvatures

at the interior points X2 through XN–1, we apply the remaining unused conditions

mentioned in (2). That is, matching the slopes of two adjacent cubic curves at these

© 2001 by CRC Press LLC

interior points. To match the slope at Xi, first we need to have the slope equation

for the preceding interval, that is from Xi–1 to Xi, for which the increment is denoted

as Hi–1. From Equation 16, we can easily write out that slope equation as:

2

?H

? ( X ? X )2 H ? Y ? Y

X i ? X) ?

(

i ?1

1

i ?1

i

?

? + Yi?? ?

Y? = Yi???1 ?

?

? i?+ i

H i ?1

2 H i ?1 ?

2 H i ?1

6?

?? 6

?

?

?

?

(17)

Using Equations 16 and 17 and matching the slopes at the interior point Xi and

after collecting terms, we obtain:

? Y ? Yi Yi ? Yi ?1 ?

H i ?1Yi???1 + 2(H i ?1 + H i )Yi?? + H i Yi??+1 = 6? i ?1

?

H i ??

? H i ?1

(18)

This equation can be applied for all interior points, that is, at X = Xi for i =

2,3,…,N–1. Hence, we have N–2 equations for solving the N–2 unknown curvatures,

(Y?)i for i = 2,3,…,N–1 when the X and Y coordinates of N + 1 points are specified

for a cubic-spline curve fit. Upon substituting the calculated curvatures into Equation

15, we obtain the desired cubic polynomial for each interval of X.

If the N given points, (Xi,Yi) for i = 1,2,…,N, has a periodic pattern for every

increment of XN-X1, we can change the above formation for the open case to suit

this particular need by requiring that the points be specified with YN = Y1 and that

curvatures also should be continuous at the first and last points. That is to remove

the 4th rule, and also one condition each for the 2nd and 3rd rules described in (2).

Equation 18 is to be used for i = 1,2,…,N to obtain N equations for solving the N

curvatures. For obtaining the first and last equations, we utilize the fact that since

Y and its derivatives are periodic, in addition to YN = Y1, (Y?)N = (Y?)1, HN = H1,

we also have YN + 1 = Y2, Y0 = YN–1, (Y?)N + 1 = (Y?)2, (Y?)0 = (Y?)N–1, and H0 = HN–1.

A program called CubeSpln has been prepared to handle both the nonperiodic

and periodic cases. It formulates the matrix equation [A]{Y?} = {C} for solving the

curvatures at Xi for i = 1,2,…,N based on Equation 18. A Gaussian elimination

scheme is needed by this subroutine for obtaining the solutions of Y?. Program

CubeSpln also has a block of statements for plotting of the spline curves. This

subroutine is listed below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Application

When program CubeSpln is run, it gives a plot of the cubic spline curves as

shown in Figure 5. The given points are marked with + symbols. Between every two

adjacent points, a third-order polynomial is derived. There are two different thirdorder polynomials for the left and right sides of every in-between points, at which

the slopes and curvatures determined by the two third-order polynomials are both

continuous.

MATLAB APPLICATION

MATLAB has a function file called spline.m which can be applied to perform

the curve fit of cubic spline. The function has three arguments: the first and second

arguments are for the coordinates of the data points for which a cubic spline curve

fit is to be obtained, and the third argument should be an array containing a more

finely spaced abscissa at which the curve ordinates should be calculated for plotting

the spline curve. Let us redo the problem for which Figure 5 has been obtained. The

MATLAB application and the resulting display are as follows:

© 2001 by CRC Press LLC

FIGURE 5. When program CubeSpln is run, it gives a plot of the cubic spline curves.

The plot of the spline curve using XC and YC data superimposed with the given

5 (X,Y) points marked by the * character is shown in Figure 6 which is identical to

Figure 5 except different in the ranges of the axes.

For dealing with data sets, more features of plot.m of MATLAB can be utilized.

Figure 7 shows how different data sets can be marked with different characters, axes

can be labeled, and various text can be added. The interactive MATLAB commands

entered are:

Notice the commands xlabel, ylabel, and text are adding labels for the x and y

axes, and text, respectively. The specific content string of label or text is to be spelled

out inside the two single quotation signs. The first two arguments for text are the

coordinates, at which the left lower corner of the first character of that string.

© 2001 by CRC Press LLC

FIGURE 6. The plot of the spline curve using XC and YC data superimposed with the given

5 (X,Y) points marked by the * character is identical to Figure 5 except different in the ranges

of the axes.

FIGURE 7. How different data sets can be marked with different characters, axes can be

labeled, and various text can be added.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

Mathematica has an interpolation function which fits a collection of specified

points by a polynomial. The command InterpolationOrder specifies the order of

polynomial and the default order is 3. Here are some examples of applications and

plots.

Input[1]: =

x = {1,2,3,4,5}

Output[1] =

{1, 2, 3, 4, 5}

Input[2]: =

y = {2,4,7,8,11}

Output[2] =

{2, 4, 7, 8, 11}

Input[3]: =

Plot[Evaluate[Interpolation[y]][x], {x,1,5},

Frame->True, AspectRatio->1]

Output[3] =

FIGURE 8.

© 2001 by CRC Press LLC

Evaluate calculates the values of the fitted polynomial within the specified

interval. In this case, the y values are interpolated using the cubic polynomial.

To add more features to a plot, Mathematica has commands Text allowing a

string of characters to be placed at a desired location, and FrameLabel allowing

labels to be added to the axes. In fact, Figure 8 is not the result of Input[3] but the

addition of the five + markers with the statement

Input[4]: =

Show[%,Graphics[Table[Text[“+”,{x[[i]],y[[i]]}],{i,1,5}]]]

Output[4] =

—Graphics—

% refers to the previous output and %% refers to the next-to-the-last output,

and so on. Show and Graphics enable additional graphics to be shown. Table lists

a series of entries specified by the running index variable which in this case is I

having a starting value of 1 and incremented by 1 (omitted), and ending at a value

equal to 5. Notice that the four * markers are not exactly located at the coordinates

specified, for the reason that the character itself is not centered but offsets like a

superscript.

As another example, consider the plot shown in Figure 9 which is resulted from

of the following statements:

Input[1]: =

X = {1,2,3,4}

Output[1] =

{1, 2, 3, 4}

Input[2]: =

Y = {2, 4, 7, 13}

Output[2] =

{2, 4, 7, 13}

Input[3]: =

X1 = {0.5, 1.2, 2.5, 3.7}

Output[3] =

{0.5, 1.2, 2.5, 3.7}

© 2001 by CRC Press LLC

Input[4]: =

Y1 = {3, 6, 5, 11}

Output[4] =

{3, 6, 5, 11}

Input[5]: =

X2 = {3.0, 3.6, 4.2, 5.1}

Output[5] =

{3.0, 3.6, 4.2, 5.1}

Input[6]: =

Y2 = {3, 6, 8, 11}

Output[6] =

{3, 6, 8, 11}

Input[7]: =

g1 = Show[Graphics[Table[Text[“*”,{X[[i]],Y[[i]] }],{i,1,5}],

Table[Text[“ + ”,{X1[[i]],Y1[[i]]}],{i,1,5}],

Table[Text[“.”,{X2[[i]],Y2[[i]]}],{i,1,5}]]]

Output[7] =

—Graphics—

Input[8]: =

g2 = Show[g1, Frame->True, AspectRatio->1,

FrameLabel->{“X”,“Y”}]

Output[8] =

—Graphics—

Input[9]: =

Show[g2,Graphics[Text[“X–Y — *”,{0.7,12},{–1,0}],

Text[“X1–Y1 — + ”,{0.7,11},{–1,0}],

Text[“X2–Y2 — .”,{0.7,10}],{–1,0}]]]

© 2001 by CRC Press LLC

Output[9] =

FIGURE 9.

The two intermediate plots designated as g1 and g2 are actually displayed on

screen but not presented here. Only the final plot showing all of the ingredients is

presented in Figure 9. Giving plot names facilitates the later referral; a better arrangement than using the % option. In Figure 9, it is also illustrated that the label for the

vertical axis is rotated by 90 degrees.

Mathematica has a package called SplineFit can also be called into service for

the need of spline curve-fit. For creating Figure 9, we may enter the request as

follows:

Input[1]: = 1]

© 2001 by CRC Press LLC

Output[4] = —Graphics—

Input[5]: = Show[%,Graphics[Table[Text[“X”,{XYS[[i]]}],{i,1,5}]]]

Output[5] = —Graphics—

In Input[1],

Exact, Least-Squares, and

Cubic Spline Curve-Fits

2.1 INTRODUCTION

Engineers conduct experiments and collect data in the laboratories. To make use of

the collected data, these data often need to be fitted with some particularly selected

curves. For example, one may want to find a parabolic equation y = c1 + c2x + c3x2

which 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 straight

line, but the reason that they are not is because of bad calibration of the measuring

equipment 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 + c2x

for the stress and strain data collected for a stretching test of a metal bar in the

elastic range, then the question of how to determine the two coefficients c1 and c2

is a matter of deciding on which criterion to adopt. The Least-Squares method is

one of the criteria which is most popularly used. The two cases cited are the

consideration 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, the

curve 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 four

particular functions, namely, fi(t) = sin(2i–1)t for i = 1,2,3,4., and to determine the

coefficients 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 can

be found to convert the equation describing the curve into linear equations. This is

discussed in a section devoted to transformed least-squares curve-fit.

Another commonly applied curve-fit technique is the cubic spline method which

allows smooth cubic equations to be derived to ensure continuous slopes and curvatures passing all given points. The mathematics involved in this method will be

presented.

In the following sections, we shall discuss the development of the programs

ExactFit, LeastSq1, LeastSqG, and CubeSpln for the four curve-fit needs mentioned above.

2.2 EXACT CURVE FIT

As another example of solving a matrix equation, let us consider the problem

of 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 of

the three points into the parabolic equation, we can obtain:

c1 + c2 x i + c3x 2i = y i

for 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 coefficient

matrix whose elements if denoted as ai,j are to be calculated using the formula:

a i, j = x ij?1

for i, j = 1, 2, 3

(3)

It is easy to extend the above argument for the general case of exactly fitting N

given points by a N-1st degree polynomial y = c1 + c2x + ••• + cNxN–1. The only

modification needed is to allow the indices i and j in Equations 1 and 3 to be extended

to reach a value of N. A program called ExactFit has been prepared for this need

by utilizing the subroutine GauJor to solve for the vector {C} from Equation 1 for

the general case of exactly fitting N given points. Listings for both FORTRAN and

QuickBASIC versions along with sample numerical results are presented below.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

© 2001 by CRC Press LLC

QUICKBASIC VERSION

Sample Application

MATLAB Application

For 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 LLC

Notice that the coefficient {C} for the curve-fit polynomial is obtained by

solving [A]{C} = {Y} where matrix [A] is formed by substituting the X values

into 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] and

MATLAB matrix operation of C = A\Y which premultiplies {Y} by [A]–1 to

obtain {C}.

Also, this exact curve-fit problem can be treated as a special case of fitting N

given points by a linear combination of N selected functions which in this case

happens to be the polynomial terms of x0 to xN–1, by the least-squares method. A m

file called LeastSqG.m which is discussed in the program LeastSqG can be readily

applied to treat such a exact curve-fit problem. Here, we demonstrate how LeastSqG.m is used by MATLAB interactive operation in solving the sample problems

previously presented in the FORTRAN and QuickBASIC versions of the program

ExactFit. First, a function must be prepared to describe the ith selected function xi

as follows:

The results of four MATLAB applications are:

© 2001 by CRC Press LLC

Notice 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. The

solution of [A]{C} = {R} is simply implemented by MATLAB with the statement

C = 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 listed

above:

Notice that for application of polyval.m, MATLAB needs the coefficients of the

polynomial arranged in descending order. Since the array C contains the coefficients

© 2001 by CRC Press LLC

FIGURE 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 the

curve values for 1?X?3 and with an increment of 0.1. Figure 1 shows that the

parabolic curve passes through all of the three given points.

2.3 PROGRAM LEASTSQ1 — LEAST-SQUARES LINEAR CURVE-FIT

Program LeastSq1 is designed for curve-fitting of a set of given data by a linear

equation based on the least-squares criterion.2 If only two points are specified, a

linear equation which geometrically describes a straight line can be uniquely derived

because the line must pass the two specified points. This is the case of exact fit. (See

programs Gauss and LagrangI for examples of exact fit.) However, the specified

data are often recorded from a certain experiment and due to inaccurate calibration

of 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 by

a type of predetermined equation. Under such a circumstance, a so-called forced fit

is then required. As a simple example, supposing that we expect the measured set

of three data points (Xi,Yi) for i = 1,2,3 to satisfy a linear law Y = c1 + c2X. If these

three points happen to fall on a straight line, then we have a case of exact fit and

© 2001 by CRC Press LLC

the coefficients c1 and c2 can be uniquely computed. If these three points are not all

on 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 selected

criterion, 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 minimize

the errors between the straight line and the points to be as small as possible. Since

the given points may be above, on, or below the straight line, and these errors may

cancel from each other if they are added together, we consider the sum of the squares

of these errors. Let us denote yi to be the value of Y at Xi and S be the sum of the

errors denoted as ei for i = 1 to N, then we can write:

N

S = e12 + e 22 + … + e 2N =

?e

2

i

(1)

i =1

where for i = 1,2,…,N

e 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 a

function of ci and c2. To find a particular set of values of c1 and c2 such that S reaches

a 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 ?

?

1

1

N

?e i

? e ?c

i

1

i =1

and

? ?e

?e ?

?e

?S

= 2? e1 1 + e 2 2 + … + e N N ? = 2

?

c

?

c

?c2

?c2 ?

?

2

2

N

?e i

? e ?c

i

i =1

2

From Equation 2, we note that ?ei/?c1 = 1 and ?ei/?c2 = Xi. Consequently, the

two extremum conditions lead to two algebraic equations for solving the coefficients

c1 and c2:

?

?

?

© 2001 by CRC Press LLC

N

?

i =1

?

?

1? c1 + ?

?

?

N

?

i =1

?

X i ? c2 =

?

N

?Y

i

i =1

(3)

and

?

?

?

N

?

i =1

?

?

X i ? c1 + ?

?

?

N

?

i =1

?

X 2i ? c2 =

?

N

?X Y

i

i

(4)

i =1

Program LeastSq1 provided in both FORTRAN and QuickBASIC versions is

developed using the above two equations. It can be readily applied for calculating

the coefficients c1 and c2. Two versions are listed and sample applications are

presented below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

MATLAB APPLICATION

A 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-squares

criterion. The function polyfit has three arguments, the first and second arguments

are the X and Y coordinate arrays of the given points, and the third argument specifies

to what degree the fitted polynomial is required. For linear fit, the third argument

should be set equal to 1. The following shows how the results obtained for the sample

problem 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.4740

If the third argument for the function polyfit is changed (from 1) to 2, 3, and

4, we also can obtain the least-squares fits of the five given points with a quadratic,

cubic, and quartic polynomials, respectively. When the third argument is set equal

to 4, we then have the case of exact curve-fit of five points by a fourth-order

polynomial. Readers are referred to the program ExactFit for more discussions.

Also, it is of interest to know whether one may select an arbitrary set of

functions and linearly combine them for least-squares fit, instead of the unbroken

set of polynomial terms X0, X1, X2, …, XN. Program LeastSqG to be presented

in the next section will discuss such generalized least-squares curve-fit. But before

we do that, let us first look into a situation where program LeastLQ1 can be

applied for a given set of data after some mathematical transformations are

employed to modify the data.

Transformed Least-Squares Curve-Fit

There are occasions when we know in advance that a given set of data supposed

to fall on a curve described by exponential equations of the type:

© 2001 by CRC Press LLC

y = b1e b2x

(5)

y = c1xe c2x

(6)

or

To determine the coefficients b1 and b2, or, c1 and c2 based on the least-squares

criterion, Equation 5 or 6 need to be first transformed into a linear form. To do so,

let us first consider Equation 5 and take natural logarithm of both sides. It gives:

n y = n b1 + n e b2x = n b1 + b 2 x

(7)

If we introduce new variable z, and new coefficients a1 and a2 such that:

z = n y,

a1 = n b1,

and

a 2 = b2

(8,9,10)

Then Equation 7 becomes:

z = a1 + a 2 x

(11)

Hence, if we need to determine the coefficients b1 and b2 for Equation 2.8 based

on N pairs of (xi,yi), for i = 1 to N, values and on the least-squares criterion, we

simply generate N z values according to Equation 2.11 and then use the N (xi,zi)

values as input for the program LeastSq1 and expect the program to calculate a1 and

a2. Equations 9 and 10 suggest that b2 is to have the value of a2 while b1 should be

equal to e raised to the a1 power, or, EXP(a1) with EXP being the exponential function

available in the FORTRAN or QuickBASIC library.

Equation 6 can be treated in a similar manner by taking logarithms of both sides

to obtain:

n y = n c1 + n x + n e c2x = n c1 + n x + c2 x

or

n y? n x = n

y

= n c1 + c2 x

x

(12)

If we introduce new variable z, and new coefficients a1 and a2 such that:

z = n(y x),

a1 = n c1,

and

a 2 = c2

(13,14,15)

Then, Equation 12 becomes Equation 11 and a1 and a2 can be obtained by the

program LeastSq1 using the data set of (xi,yi/xi) values.

As a numerical example, consider the case of a set of nine stress-versus-strain

( vs. ) data collected from a stretching test of a long bar: (.265,1025), (.4,1400),

(.5,1710), (.7,2080), (.95,2425), (1.36,2760), (2.08,3005), (2.45,2850), and

© 2001 by CRC Press LLC

(2.94,2675) where the units for the strains and stresses are in microinch x 102 and

lb/in2, respectively. When program LeastSq1 is applied for the modified data of

(,), the coefficients for the linear fit are a1 = 15.288 and a2 = –537.71. Consequently, according to Equations 13 and 14, and realizing that x and y are now and

, respectively, we have arrived at = 4.3615 ? 106e–537.71. The derived curve and

the given points are plotted in Figure 2 which shows the curve passes the origin as

it should.

FIGURE 2. The curve passes the origin as it should.

© 2001 by CRC Press LLC

2.4 PROGRAM LEASTSQG — GENERALIZED

LEAST-SQUARES CURVEFIT

Program LeastSqG is designed for curve-fitting of a set of given data by a linear

combination of a set of selected functions based on the least-squares criterion.2

Let us consider N points whose coordinates are (Xk,Yk) for k = 1 to N and let

the M selected function be f1(X) to fM(X) and the equation determined by the leastsquares curve-fit be:

M

Y( X) = a1f1 ( X) + a 2 f2 ( X) + … + a M fM ( X) =

? a f (X)

j j

(1)

j=1

The least-squares curve-fit for finding the coefficients c1–M is to minimize the

errors between the computed Y values based on Equation 1 and the given Y values

at all Xk’s for k = 1 to N. Let us denote yk to be the value of Y calculated at Xk

using Equation 1 and S be the sum of the errors denoted as ek for k = 1 to N. Since

the yk could either be greater or less than Yk, these errors ek’s may cancel from each

other if they are added together. We therefore consider the sum of the squares of

these errors and write:

N

S = e12 + e 22 + … + e 2N =

?e

2

k

(2)

k =1

where for k = 1,2,…,N

?

e k = y k ? Yk = ?

?

?

M

?

j=1

?

a jfj ( X k )? ? Yk

?

?

(3)

It is obvious that since Xk and Yk are constants, the sum S of the errors is a

function of a1 to M. To find a particular set of values of a1 to M such that S reaches a

minimum, we may therefore base on calculus3 and utilize the conditions ?S/?ai = 0

for i = 1 to M. From Equation 2, we can have the partial derivatives of S with respect

to ai’s, for i = 1 to M, as:

N

? ?e

?e ?

?e

?e

?S

= 2? e1 1 + e 2 2 + … + e N N ? = 2

ek k

?a i

?a i

?a i

?a i ?

? ?a1

k =1

?

From Equation 3, we note that ?ek/?ai = fi(Xk). Consequently, the M extremum

conditions, ?S/?ai = 0 for i = 1 to M, lead to M algebraic equations for solving the

coefficients a1 to M. That is, for i = 1 to M:

M

?

? N

fi ( X k )fj ( X k )?a j =

?

??

?? k =1

??

j=1

© 2001 by CRC Press LLC

N

? f (X )Y

i

k =1

k

k

(4)

If we express Equation 4 in matrix notation, it has the simple form:

[C]{A} = {R}

(5)

where [C] is a MxM square coefficient matrix, and {A} and {R} are Mx1 column

matrices (vectors). {A} contains the unknown coefficients a1 to M needed in Equation 1.

If we denote the elements in [C] and {R} as cij and ri, respectively, Equation 5

indicates that these elements are to be calculated using the formulas, for i = 1,2,…,M:

N

cij =

? f (X )f (X )

i

k

j

k

(6)

k =1

and

N

ri =

? f (X )Y

i

k

k

(7)

k =1

The above derivation appears to be too mathematical; a few examples of actual

curve-fit will clarify the procedure involved. As a first example, consider the case

of selecting two (M = 2) functions f1(X) = 1 and f2(X) = X to fit three given points

(N = 3), (X1,Y1) = (1,1), (X2,Y2) = (2.6,2), and (X3,Y3) = (2.8,2). Equations 6 and 7

then provide the following:

and

Hence, the system of two linear algebraic equations for finding a1 and a2 for the

straight-line equation is:

?3

?6.4

?

© 2001 by CRC Press LLC

6.4 ? ? a1 ? ? 5 ?

=

15.6?? ??a 2 ?? ??11.8??

The solution can be obtained by application of Cramer’s Rule to be a1 = 0.42466

and a2 = 0.58219. More examples will be given after we discuss how computer

programs can be written to compute [C] and {R} and then solve for {A}.

Program LeastSqG provided in both FORTRAN and QuickBASIC versions

is developed using the above two equations. It can be readily applied for calculating

the coefficients c1 to M. Both QucikBASIC and FORTRAN versions are listed and

sample applications are presented below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Applications

When four functions are selected as those listed in SUB FS, an interactive

application of the program LeastSqG QuickBASIC version using the input data

entered through keyboard has resulted in a screen display of:

If three sinusoidal functions, sin(x/20), sin(3x/20), and sin(5x/20) were

selected and replacing those listed in SUB FS, another interactive application of the

program LeastSqG QuickBASIC version is shown below.

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

By selecting the four functions listed in Subroutine FS, an interactive application

of the program LeastSqG using the input data given below has resulted in a screen

display of:

MATLAB APPLICATION

A LeastSqG.m file can be created and added to MATLAB m files which will

take N sets of X and Y points and fitted by a linear combination of M selected

functions in the least-squares sense. The selected functions can be specified in

another m file called FS.m (using the same name as in the FORTRAN and QuickBASIC versions). First, let us look at a version of LeastSqG.m:

Notice that the coefficients C’s is obtained by solving [A]{C} = {R} as in the

text. For MATLAB, a simple matrix multiplication of the inverse of [A] to and on

the left of the vector {R}. Complete execution of LeastSqG.m will be indicated by

a display of the matrix [A], vector {R}, and the solution vector {C}. The expression

feval(funfcn,i,X(k)) in the above program is to evaluate the ith function at X(k)

defined in a function file to be specified when LeastSqG.m is applied which is to

be illustrated next.

© 2001 by CRC Press LLC

Consider the case of given 5 (X,Y) points (N = 5) which are (1.4,2.25), (3.2,15),

(4.8,26.25), (8,33), and (10,35). And, the selected functions are sin(X/20),

sin(3X/20), and sin(5X/20). That is, M = 3. The supporting function FS.m is

simply:

function F=FS(i,xv)

F=sin((2.*i-1).*xv.(pi./20);

Having prepared this file FS.m on a disk which is in drive A, we can now apply

LeastSqG.m interactively as follows:

>> X = [1.4,3.2,4.8,8,10]; Y = [2.25,15,26.25,33,35]; C = LeastSqG(A:FS,X,Y,3,5)

The results shown on screen are:

If four functions X, X2, sin(X), and cos(X) are selected, we may change the

FS.m file to:

© 2001 by CRC Press LLC

Same as for the FORTRAN and QuickBASIC versions, if we are given six

(X,Y) points, (1,2), (2,4), (3,7), (4,11), (5,23), and (6,45), MATLAB application of

LeastSqG.m will be:

>> X=[1,2,3,4,5,6]; Y=[2,4,7,11,23,45]; C-LeastSqG('A:FS',X,Y,4,6)

The results are:

Notice that the values in [A] and {R} are to be multiplied by the factor 1.0e +

003 as indicated. For saving space, [A], {R}, and {C} are listed side-by-side when

actually they are displayed from top-to-bottom on screen. To further utilize the

capability of MATLAB, the obtained {C} values are checked to plot the fitted curve

against the provided six (X,Y) points. The following additional statements are

entered in order to have a screen display as illustrated in Figure 3:

FIGURE 3. The ‘*’ argument in the second plot statement requests that the character * is

to be used for plotting the given points,

© 2001 by CRC Press LLC

The statement XC = [1:0.2:6] generates a vector of XC containing value from

1 to 6 with an increment of 0.2. The command “hold” enables the first plot of XC

vs. YC which is the resulting least-squares fitted curve, to be held on screen until

the given points (X,Y) are superimposed. The ‘*’ argument in the second plot

statement requests that the character * is to be used for plotting the given points, as

illustrated in Figure 3.

Next, another example is given to show how MATLAB statements can be applied

directly with defining a m function, such as FS which describes the selected set of

functions for least-squares curve fit. Consider the problem of least-squares fit of the

points (1,2), (3,5), and (4,13) by the linear combination Y = C1f1(X) + C2f2(X) where

f1(X) = X–1 and f2(X) = X2. An interactive application of MATLAB may go as follows:

The resulting curve is plotted in Figure 4 using 31 points, (XC,YC), calculated

based on the equation Y = –7.0526(X–1) + 2.1316X2 for X values ranging from 1

to 4. The three given points, (X,Y), are superimposed on the graph using ‘*’

character. Notice from the above statements, the coefficients C(1) and C(2) are solved

from the matrix equation [A]{C} = {R} where the elements in [A] are generated

using interactively entered statement and so are the elements of {R}. MATLAB

matrix operations such as transposition, subtraction, multiplication, and inversion

are all involved. Also, notice that here no use of MATLAB ‘hold’ as for generating

Figure 1, is necessary if X, Y, XC, and YC are all used as arguments in calling the

© 2001 by CRC Press LLC

FIGURE 4. The curve is plotted using 31 points, (XC,YC), calculated based on the equation

Y = –7.0526(X–1) + 2.1316X2 for X values ranging from 1 to 4.

plot function. The statement XC = 1:0.1:4 generates a XC row matrix having elements starting from a value equal to 1, equally incremented by 0.1 until a value of

4 is reached.

MATHEMATICA APPLICATIONS

Mathematica has a function called Fit which least-squares fits a given set of

(x,y) points by a linear combination of a number of selected functions. As the first

example of interactive application, let us find the best straight line which gives the

least squared errors for fitting a set of five points. This is the case of two selected

functions f1(x) = 1 and f2(x) = x. the interactive application goes as follows:

In[1]: = Fit[{{1,2},{2,5},{3,8},{5,11},{8,24}}, {1, x}, x]

Out[1]: = –1.47403 + 3.01948 x

Notice that Fit has three arguments: first argument specifies the data set, second

argument lists the selected function, and the third argument is the variable for the

© 2001 by CRC Press LLC

derived equation. In case that three points are given and the two selected functions

are f1(x) = x–1 and f2(x) = x2, then the interactive operation goes as follows:

In[2]: = Fit[{{1,2},{3,5},{4,13}}, {x–1, x^2}, x]

Out[2]: = –7.05263 (–1 + x) + 2.13158 x

Two other examples previously presented in the MATLAB applications can also

be considered and the results are:

In[3]: = (Fit[{{1,2}, {2,4}, {3,7}, {4,11}, {5,23}, {6,45}},

{x, x^2, Sin[x], Cos[x]}, x])

Out[3]: = –4.75756 x + 2.11159 x2 – 0.986915 Cos[x] + 5.76573 Sin[x]

and

In[4]: = (Fit[{{1.4, 2.25}, {3.2, 15}, {4.8, 26.25}, {8, 33}, {10, 35}},

{Sin[Pi*x/20], Sin[3*Pi*x/20], Sin[5*Pi*x/20]}, x])

Pi x ?

3 Pi x ?

5 Pi x ?

Out[4]: = 35.9251 Sin ??

? 1.19261 Sin ??

? 3.46705 Sin ??

? 20 ??

? 20 ??

? 20 ??

All of the results obtained here are in agreement with those presented earlier.

2.5 PROGRAM CUBESPLN — CURVE FITTING

WITH CUBIC SPLINE

If a set of N given points (Xi,Yi) for i = 1, 2,…,N is to be fitted with a polynomial

of N–1 degree passing these points exactly, the polynomial curve will have many

fluctuations between data points as the number of points, N, increases. A popular

method for avoiding such over-fluctuation is to fit every two adjacent points with a

cubic equation and to impose continuity conditions at the point shared by two

neighboring cubic equations. This is like using a draftsman’s spline attempting to

pass through all of the N given points. It is known as cubic-spline curve fit. For a

typical interval of X, say Xi to Xi + 1, the cubic equation can be written as:

Y = a i X3 + bi X 2 + ci X + d i

(1)

where ai, bi, ci, and di are unknown coefficients. If there are N–1 intervals and each

interval is fitted by a cubic equation, the total number of unknown coefficients is

equal to 4(N–1). They are to be determined by the following conditions:

1. Each curve passes two specified points. The first and last specified points

are used once by the first and N-1st curves, respectively, whereas all

interior points are used twice by the two cubic curves meeting there. This

gives 2 + 2(N–2), or, 2N–2 conditions.

© 2001 by CRC Press LLC

2. Every two adjacent cubic curves should have equal slope (first derivative

with respect to X) at the interior point which they share. This gives N–2

conditions.

(2)

3. For further smoothness of the curve fit, every two adjacent cubic curves

should have equal curvature (second derivative with respect to X) at the

interior point which they share. This gives N–2 conditions.

4. At the end points, X1 and XN, the nature spline requires that the curvatures

be equal to zero. This gives 2 conditions.

Instead of solving the coefficients in Equation 1, the usual approach is to apply

the above-listed conditions for finding the curvatures, Y? at the interior points, that

is for i = 2,3,…,N–1 since (Y?)1 = (Y?)N = 0. To do so, we notice that if Y is a cubic

polynomial of X, Y? is then linear in X and can be expressed as:

Y?? = AX + B

(3)

If this is used to fit the ith interval, for which the increment in X is here denoted

as Hi = Xi + 1–Xi, we may replace the unknown coefficients A and B with the

curvatures at Xi and Xi + 1, (Y?)i and (Y?)i + 1 by solving the two equations:

Yi?? = AX i + B

and

Yi??+1 = AX i +1 + B

(4,5)

By using the Cramer’s rule, it is easy to obtain:

A=

Yi?+1 ? Yi??

X i +1 ? X i

and

B = Y?? ?

X i (Yi??+1 ? Yi??)

X i +1 ? X i

(6,7)

Consequently, Equation 3 can be written as:

Y?? =

Yi??

(X ? X) + YHi??+1 (X ? Xi )

H i i +1

i

(8)

Equation 8 can be integrated successively to obtain the expressions for Y and

Y as:

2

2

? Yi??

Y??

Y? +

X i +1 ? X) + i +1 ( X ? X i ) + C1

(9)

(

2Hi

2H i

and

3

3

Y??

Y??

Y = i ( X i +1 ? X) + i +1 ( X ? X i ) + C1X + C2

(10)

6H i

6H i

© 2001 by CRC Press LLC

The integration constants C1 and C2 can be determined by the conditions that at

Xi, Y = Yi and at Xi + 1, Y = Yi + 1. Based on Equation 10, the two conditions lead to:

Yi =

Yi?? 2

H + C1X i + C2

6 i

and

Yi +1 =

Yi??+1 2

H + C1X i +1 + C2

6 i

(11,12)

Again, Cramer’s rule can be applied to obtain:

C1 =

(Y??? Y?? )H

i +1

i

6

i

?

Yi ? Yi +1

Hi

(13)

and

C2 = ?

(Xi+1Yi??? XiYi??+1 )Hi ? XiYi+1 ? Xi+1Yi

6

Hi

Hi

(14)

Substituting C1 and C2 into Equations 11 and 12 and rearranging into terms

involving the unknown curvatures, the expressions for the cubic curve are:

Y=

3

? Y?? ? ( X ? X )3

?

Yi?? ? ( X i +1 ? X)

i

?

? H i ( X i +1 ? X)? + i +1 ?

? H i ( X ? X i )?

Hi

6 ?

6 ? Hi

??

??

?

?

? X ? X?

? X ? Xi ?

+ Yi ? i +1

+ Yi +1 ?

?

?

? Hi ?

? Hi ?

(15)

? H ( X ? X )2 ?

? ( X ? X )2 H ? Y ? Y

i

i

? + Yi??+1 ?

Y? = Yi?? ? i ? i +1

? i ? + i +1

Hi

2Hi

6?

??

?? 6

?? 2 H i

?

(16)

and

Equation 15 indicates that the cubic curve for the ith interval is completely

defined if in addition to the specified values of Yi and Yi + 1, the curvatures at Xi and

Xi + 1, (Y?)i and (Y?)i + 1 respectively, also can be found. To find all of the curvatures

at the interior points X2 through XN–1, we apply the remaining unused conditions

mentioned in (2). That is, matching the slopes of two adjacent cubic curves at these

© 2001 by CRC Press LLC

interior points. To match the slope at Xi, first we need to have the slope equation

for the preceding interval, that is from Xi–1 to Xi, for which the increment is denoted

as Hi–1. From Equation 16, we can easily write out that slope equation as:

2

?H

? ( X ? X )2 H ? Y ? Y

X i ? X) ?

(

i ?1

1

i ?1

i

?

? + Yi?? ?

Y? = Yi???1 ?

?

? i?+ i

H i ?1

2 H i ?1 ?

2 H i ?1

6?

?? 6

?

?

?

?

(17)

Using Equations 16 and 17 and matching the slopes at the interior point Xi and

after collecting terms, we obtain:

? Y ? Yi Yi ? Yi ?1 ?

H i ?1Yi???1 + 2(H i ?1 + H i )Yi?? + H i Yi??+1 = 6? i ?1

?

H i ??

? H i ?1

(18)

This equation can be applied for all interior points, that is, at X = Xi for i =

2,3,…,N–1. Hence, we have N–2 equations for solving the N–2 unknown curvatures,

(Y?)i for i = 2,3,…,N–1 when the X and Y coordinates of N + 1 points are specified

for a cubic-spline curve fit. Upon substituting the calculated curvatures into Equation

15, we obtain the desired cubic polynomial for each interval of X.

If the N given points, (Xi,Yi) for i = 1,2,…,N, has a periodic pattern for every

increment of XN-X1, we can change the above formation for the open case to suit

this particular need by requiring that the points be specified with YN = Y1 and that

curvatures also should be continuous at the first and last points. That is to remove

the 4th rule, and also one condition each for the 2nd and 3rd rules described in (2).

Equation 18 is to be used for i = 1,2,…,N to obtain N equations for solving the N

curvatures. For obtaining the first and last equations, we utilize the fact that since

Y and its derivatives are periodic, in addition to YN = Y1, (Y?)N = (Y?)1, HN = H1,

we also have YN + 1 = Y2, Y0 = YN–1, (Y?)N + 1 = (Y?)2, (Y?)0 = (Y?)N–1, and H0 = HN–1.

A program called CubeSpln has been prepared to handle both the nonperiodic

and periodic cases. It formulates the matrix equation [A]{Y?} = {C} for solving the

curvatures at Xi for i = 1,2,…,N based on Equation 18. A Gaussian elimination

scheme is needed by this subroutine for obtaining the solutions of Y?. Program

CubeSpln also has a block of statements for plotting of the spline curves. This

subroutine is listed below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Application

When program CubeSpln is run, it gives a plot of the cubic spline curves as

shown in Figure 5. The given points are marked with + symbols. Between every two

adjacent points, a third-order polynomial is derived. There are two different thirdorder polynomials for the left and right sides of every in-between points, at which

the slopes and curvatures determined by the two third-order polynomials are both

continuous.

MATLAB APPLICATION

MATLAB has a function file called spline.m which can be applied to perform

the curve fit of cubic spline. The function has three arguments: the first and second

arguments are for the coordinates of the data points for which a cubic spline curve

fit is to be obtained, and the third argument should be an array containing a more

finely spaced abscissa at which the curve ordinates should be calculated for plotting

the spline curve. Let us redo the problem for which Figure 5 has been obtained. The

MATLAB application and the resulting display are as follows:

© 2001 by CRC Press LLC

FIGURE 5. When program CubeSpln is run, it gives a plot of the cubic spline curves.

The plot of the spline curve using XC and YC data superimposed with the given

5 (X,Y) points marked by the * character is shown in Figure 6 which is identical to

Figure 5 except different in the ranges of the axes.

For dealing with data sets, more features of plot.m of MATLAB can be utilized.

Figure 7 shows how different data sets can be marked with different characters, axes

can be labeled, and various text can be added. The interactive MATLAB commands

entered are:

Notice the commands xlabel, ylabel, and text are adding labels for the x and y

axes, and text, respectively. The specific content string of label or text is to be spelled

out inside the two single quotation signs. The first two arguments for text are the

coordinates, at which the left lower corner of the first character of that string.

© 2001 by CRC Press LLC

FIGURE 6. The plot of the spline curve using XC and YC data superimposed with the given

5 (X,Y) points marked by the * character is identical to Figure 5 except different in the ranges

of the axes.

FIGURE 7. How different data sets can be marked with different characters, axes can be

labeled, and various text can be added.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

Mathematica has an interpolation function which fits a collection of specified

points by a polynomial. The command InterpolationOrder specifies the order of

polynomial and the default order is 3. Here are some examples of applications and

plots.

Input[1]: =

x = {1,2,3,4,5}

Output[1] =

{1, 2, 3, 4, 5}

Input[2]: =

y = {2,4,7,8,11}

Output[2] =

{2, 4, 7, 8, 11}

Input[3]: =

Plot[Evaluate[Interpolation[y]][x], {x,1,5},

Frame->True, AspectRatio->1]

Output[3] =

FIGURE 8.

© 2001 by CRC Press LLC

Evaluate calculates the values of the fitted polynomial within the specified

interval. In this case, the y values are interpolated using the cubic polynomial.

To add more features to a plot, Mathematica has commands Text allowing a

string of characters to be placed at a desired location, and FrameLabel allowing

labels to be added to the axes. In fact, Figure 8 is not the result of Input[3] but the

addition of the five + markers with the statement

Input[4]: =

Show[%,Graphics[Table[Text[“+”,{x[[i]],y[[i]]}],{i,1,5}]]]

Output[4] =

—Graphics—

% refers to the previous output and %% refers to the next-to-the-last output,

and so on. Show and Graphics enable additional graphics to be shown. Table lists

a series of entries specified by the running index variable which in this case is I

having a starting value of 1 and incremented by 1 (omitted), and ending at a value

equal to 5. Notice that the four * markers are not exactly located at the coordinates

specified, for the reason that the character itself is not centered but offsets like a

superscript.

As another example, consider the plot shown in Figure 9 which is resulted from

of the following statements:

Input[1]: =

X = {1,2,3,4}

Output[1] =

{1, 2, 3, 4}

Input[2]: =

Y = {2, 4, 7, 13}

Output[2] =

{2, 4, 7, 13}

Input[3]: =

X1 = {0.5, 1.2, 2.5, 3.7}

Output[3] =

{0.5, 1.2, 2.5, 3.7}

© 2001 by CRC Press LLC

Input[4]: =

Y1 = {3, 6, 5, 11}

Output[4] =

{3, 6, 5, 11}

Input[5]: =

X2 = {3.0, 3.6, 4.2, 5.1}

Output[5] =

{3.0, 3.6, 4.2, 5.1}

Input[6]: =

Y2 = {3, 6, 8, 11}

Output[6] =

{3, 6, 8, 11}

Input[7]: =

g1 = Show[Graphics[Table[Text[“*”,{X[[i]],Y[[i]] }],{i,1,5}],

Table[Text[“ + ”,{X1[[i]],Y1[[i]]}],{i,1,5}],

Table[Text[“.”,{X2[[i]],Y2[[i]]}],{i,1,5}]]]

Output[7] =

—Graphics—

Input[8]: =

g2 = Show[g1, Frame->True, AspectRatio->1,

FrameLabel->{“X”,“Y”}]

Output[8] =

—Graphics—

Input[9]: =

Show[g2,Graphics[Text[“X–Y — *”,{0.7,12},{–1,0}],

Text[“X1–Y1 — + ”,{0.7,11},{–1,0}],

Text[“X2–Y2 — .”,{0.7,10}],{–1,0}]]]

© 2001 by CRC Press LLC

Output[9] =

FIGURE 9.

The two intermediate plots designated as g1 and g2 are actually displayed on

screen but not presented here. Only the final plot showing all of the ingredients is

presented in Figure 9. Giving plot names facilitates the later referral; a better arrangement than using the % option. In Figure 9, it is also illustrated that the label for the

vertical axis is rotated by 90 degrees.

Mathematica has a package called SplineFit can also be called into service for

the need of spline curve-fit. For creating Figure 9, we may enter the request as

follows:

Input[1]: = 1]

© 2001 by CRC Press LLC

Output[4] = —Graphics—

Input[5]: = Show[%,Graphics[Table[Text[“X”,{XYS[[i]]}],{i,1,5}]]]

Output[5] = —Graphics—

In Input[1],

### Случайные PDF

Файл |

2.pdf |

1.pdf |

Тестирование ПО (реферат).pdf |

Теория.pdf |

Посторонний.pdf |