Azərbaycan Dövlət Neft

və Sənaye Universiteti
Numerical
Methods I
Curve-Fitting:
Lagrange Interpolation
The Lagrange interpolating polynomial is
simply a reformulation of the Newton
polynomial that avoids the computation of
divided differences.

It can be represented concisely as

For example the linear version () is

The second-order version () is

A visual depiction of the rationale behind the
Lagrange polynomial is illustrated in .

This figure shows a second-

order case.

Each of the three terms in

Eq. passes through one of
the data points and is zero at
the other two.
The summation of the three terms must,
therefore, be the unique second-order
polynomial that passes exactly through the
three points.

Equation can be very simply programmed for

implementation on a computer.

below shows pseudocode that can be employed

for this purpose.
. This algorithm is set up to compute a single th-order
prediction; is the number of data points.
The Lagrange interpolating polynomial can be
derived directly from Newtons formulation we
will do this for the first-order case only.

To derive the Lagrange form we reformulate the

divided differences.

For example the first divided difference can be

reformulated as (called the )
Substituting the last equation into

yields

Finally grouping similar terms and simplifying

yields the Lagrange form
Each term will be at and at all other sample
points.

Thus each product takes on the value of at the

sample point .

Consequently, the summation of all the

products designated by Eq. is the unique
order polynomial that passes exactly through
th

all data points.

.
Use a Lagrange interpolating polynomial of the
first and second orders to evaluate on the basis
of the following data:
.
First-order Lagrange interpolating polynomial:
.

Second-order Lagrange interpolating polynomial:

Note that, as with Newton’s method, the
Lagrange version has an estimated error of

Thus if an additional point is available at an

error estimate can be obtained.

However because the finite divided differences

are not employed as part of the Lagrange
algorithm this is rarely done.
In summary, for cases where the order of the
polynomial is unknown, the Newton method
has advantages because of the insight it
provides into the behavior of the different-order
formulas.

In addition, the error estimate represented by

Eq. can usually be integrated easily into the
Newton computation because the estimate
employs a finite difference.
Thus, for exploratory computations, Newton’s
method is often preferable.

When only one interpolation is to be performed,

the Lagrange and Newton formulations require
comparable computational effort.

However, the Lagrange version is somewhat

easier to program.
Because it does not require computation and
storage of divided differences, the Lagrange
form is often used when the order of the
polynomial is known a priori.

.
We can use the algorithm from to study a trend
analysis problem associated with our now-
familiar falling parachutist.
.
Assume that we have instrumentation to
measure the velocity of the parachutist.

The measured data obtained for a particular test

case are presented
in the table:
.
Our problem is to estimate the velocity of the
parachutist at to fill in the large gap in the
measurements between and .

We are aware that the behavior of interpolating

polynomials can be unexpected.

Therefore, we will construct polynomials of

orders and and compare the results.
. The Lagrange algorithm can be used to
construct fourth-, third-, second-, and first-order
interpolating polynomials.

The fourth-order polynomial and the input data

can be plotted
as shown in
.
. It is evident from this plot that the estimated
value of at is higher than the overall trend of
these data.

through shows plots of the results of the

computations for third-, second-, and first-order
interpolating polynomials, respectively.

It is noted that the lower the order, the lower the

estimated value of the velocity at .
. Third-order polynomial approximation
. Second-order polynomial approximation
. First-order polynomial approximation
. The plots of the interpolating polynomials indicate
that the higher-order polynomials tend to overshoot
the trend of these data.

This suggests that the first- or second-order versions

are most appropriate for this particular trend analysis.

It should be remembered, however, that because we are

dealing with uncertain data, regression would actually
be more appropriate.
The preceding example illustrates that higher-
order polynomials tend to be ill-conditioned,
that is, they tend to be highly sensitive to
round-off error.

The same problem applies to higher-order

polynomial regression.

Double-precision arithmetic sometimes helps

mitigate the problem.
However, as the order increases, there will come
a point at which round-off error will interfere
with the ability to interpolate using the simple
approaches covered to this point.
Coefficients of Interpolating
Polynomials
Although both the Newton and the Lagrange
polynomials are well-suited for determining
intermediate values between points they do not
provide a convenient polynomial of the
conventional form
A straightforward method for computing the
coefficients of this polynomial is based on the
fact that data points are required to determine
the coefficients.
Thus simultaneous linear algebraic equations
can be used to calculate the .

For example, suppose that you desired to

compute the coefficients of the parabola
Three data points are required:

Each can be substituted into Eq. to give

Thus for this case the are the knowns and the
are the unknowns.
Because there are the same number of equations
as unknowns Eq. could be solved by an
elimination method (e.g. Gaussian elimination
method).

It should be noted that the foregoing approach

is not the most efficient method that is available
to determine the coefficients of an interpolating
polynomial.
Remark
Systems such as Eq. are notoriously ill-
conditioned.

Whether they are solved with an elimination

method or with a more efficient algorithm the
resulting coefficients can be highly inaccurate
particularly for large .

When used for a subsequent interpolation they

often yield erroneous results.
In summary, if you are interested in
determining an intermediate point, employ
Newton or Lagrange interpolation.

If you must determine an equation of the form

of Eq. limit yourself to lower-order
polynomials and check your results carefully.
Inverse Interpolation
The and values in most interpolation contexts
are the dependent and independent variables
respectively.

As a consequence the values of the s are

typically uniformly spaced.

As a simple example consider a table of values

derived for the function
Now suppose that you must use the same data
but you are given a value for and must
determine the corresponding value of .
For instance for the data above suppose that you
were asked to determine the value of that
corresponded to .

For this case because the function is available

and easy to manipulate the correct answer can
be determined directly as .
Such a problem is called .

For a more complicated case you might be

tempted to switch the and values i.e. merely
plot versus and use an approach like Lagrange
interpolation to determine the result.

Unfortunately when you reverse the variables

there is no guarantee that the values along the
new abscissa the s will be evenly spaced.
For example for the result of inversion results
in the following table:

Such nonuniform spacing on the abscissa often

leads to oscillations in the resulting
interpolating polynomial this can occur even for
lower-order polynomials.
An alternative strategy is to fit an order
th

interpolating polynomial to the original data

i.e. with versus .

In most cases because the s are evenly spaced

this polynomial will not be ill-conditioned.

The answer to our problem then amounts to

finding the value of that makes this polynomial
equal to the given .
Thus the interpolation problem reduces to a
problem.

For example for the problem outlined above a

simple approach would be to fit a quadratic
polynomial to the three points:

The result would be

The answer to the inverse interpolation
problem of finding the corresponding to
would therefore involve determining the root of

For this simple case the quadratic formula can

be used to calculate
Thus, the second root is a good approximation
of the true value of .

If additional accuracy were desired a third or

fourthorder polynomial along with one of the
root location methods could be employed.
Interpolation with Equally
Spaced Data
If data is equally spaced and in ascending order
then the independent variable takes on values:

where is the or .
On this basis the finite divided differences can
be expressed in concise form.

For example, the second forward divided

difference is

which can be expressed as

That is possible because

Now recall that the is an expression of the form

A uses the function values at and , instead of

the values at and , i.e. .
Finally, the is given by

Then the is equal to

More generally, the th

order forward difference
is given by
These equations use after
the summation sign shown
as .

Each row of Pascal's triangle provides the

coefficient for each value of .
Therefore Eq. can be represented as
And in general

Using Eq. we can express Newton’s

interpolating polynomial for the case of
equispaced data as
Here the remainder is the same as

Equation is known as or the .

It can be simplified further by defining a new

quantity as follows:
This definition can be used to develop the
following simplified expressions for the terms
in Eq. :
Substituting these values into Eq. give us:

where

In addition to the forward formula and

Newton-Gregory formulas are also available.
Because both the Newton and Lagrange
polynomials are compatible with arbitrarily
spaced data you might wonder why we address
the special case of equally spaced data.

They are commonly needed to derive numerical

integration formulas that typically employ
equispaced data numerical differentiation and
integration.
And the numerical integration formulas have
relevance to the solution of ordinary differential
equations.
Extrapolation
is the process of estimating a value of that lies
outside the range of the known base points .

Note that the most accurate interpolation is

usually obtained when the unknown lies near
the center of the base points.
Obviously this is violated when the unknown
lies outside the range and consequently the
error in extrapolation can be very large.
As depicted in the the open-ended nature of
extrapolation represents a step into

the unknown
because the process
extends the curve
beyond the known
region.
As such the true curve could easily diverge from
the prediction.

Extreme care should therefore be exercised

whenever a case arises where one must
extrapolate.
.
The population in millions of the United States
of America from to can be tabulated as follows:

Fit a seventh-order polynomial to the first

points ( to ). Use it to compute the population in
by extrapolation and compare your prediction
with the actual result.
. First, the data can be entered in Matlab® as
follows:
>> t = 1920:10:1990;
>> p = [106.46 123.08 132.12 152.27
180.67 205.05 227.23 249.46]

The function can be used to compute the

coefficients as follows:
>> pol = polyfit(t,p,7);
. However, when this is done, the following
message is displayed:
Warning: Polynomial is badly conditioned.
Add points with distinct X values, reduce
the degree of the polynomial, or try
centering and scaling as described in HELP
POLYFIT.
We can follow Matlab’s suggestion by scaling
and centering the data values as in:
>> ts = (t – 1955)/35;
. Now works without an error message:
pol = polyfit(ts,p,7);

We can then use the polynomial coefficients

along with the function to predict the
population in as follows:
>> polyval(pol,(2000 - 1955)/35)
ans =
175.0800
. This is much lower than the true value of .
Insight into the problem can be gained by
generating a plot of the data and the
polynomial:
>> tt = linspace(1920,2000);
>> polp = polyval(pol,(tt -
1955)/35);
>> plot(t,p,'o', tt,polp); grid on;
xlabel('time');
. Use of a seventh-order polynomial to make a
prediction of U.S. population in based on data
from through . 250

200
population

150

100
1920 1930 1940 1950 1960 1970 1980 1990 2000
time
. As illustrated in the last , the result indicates
that the polynomial seems to fit the data nicely
from to .

However, once we move beyond the range of

the data into the realm of extrapolation, the
seventh-order polynomial plunges to the
erroneous prediction in .
Although “more is better” in many contexts, it is
absolutely not true for polynomial
interpolation.

Higher-order polynomials tend to be very ill-

conditioned — that is, they tend to be highly
sensitive to round off error.
The following example illustrates this point
nicely.
.

In , Carl Runge published a study on the

dangers of higher-order polynomial
interpolation. He looked at the following
simple-looking function:

which is now called .

. He took equidistantly spaced data points from
this function over the interval .

He then used interpolating polynomials of

increasing order and found that as he took more
points, the polynomials and the original curve
differed considerably.

Further, the situation deteriorated greatly as the

order was increased.
. Duplicate Runge’s result by using the and
functions to fit fourth- and tenth-order
polynomials to and equally spaced points
generated with Eq. .

Create plots of your results along with the

sampled values and the complete Runge’s
function.
. The five equally spaced data points can be
generated as in:
>> x = linspace(-1,1,5);
>> y = 1./(1 + 25*x.^2);

Next, a more finally spaced vector of values can

be computed so that we can create a smooth plot
of the results:
>> xx = linspace(-1,1);
.
The function can be used to generate the
coefficients of the fourth-order polynomial, and
the function can be used to generate the
polynomial interpolation at the finely spaced
values of :
>> p = polyfit(x,y,4);
>> y4 = polyval(p,xx);
. Finally, we can generate values for Runge’s
function itself and plot them along with the
polynomial fit and the sampled data:
>> yr = 1./(1 + 25*xx.^2);
>> plot(x,y,'o', xx,y4, xx,yr,'--');
grid on; xlabel('x');
ylabel('y');
. Comparison of Runge’s function (dashed line)
with a fourth-order polynomial fit to points
sampled from the function.
1

0.8

As seen in , the
0.6

0.4

polynomial does a
y

0.2

poor job of 0

following Runge’s -0.2

function. -0.4
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
x
.
Continuing with the analysis, the tenth-order
polynomial can be generated and plotted with:
>> x = linspace(-1,1,11);
>> y = 1./(1 + 25*x.^2);
>> xx = linspace(-1,1);
>> p = polyfit(x,y,10);
>> y10 = polyval(p,xx);
>> yr = 1./(1 + 25*xx.^2);
.
>> plot(x,y,'o', xx,y10,
xx,yr,'--');
Comparison of grid 2
on; xlabel('x');
Runge’s function
ylabel('y'); 1.5

(dashed line) with 1

a tenth-order
y

0.5

polynomial fit to
0

points sampled
from the function.
-0.5
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
x
As can be seen in the Figure, the fit has gotten
even worse, particularly at the endpoints!

Although there may be certain contexts where

higher-order polynomials are necessary, they
are usually to be avoided. In most engineering
and scientific contexts, lower-order polynomials
of the type described in this chapter can be used
effectively to capture the curving trends of data
without suffering from oscillations.
Thank you for attention!
text

text

text
. text
text

text
text