Lagrange’s Interpolation Polynomial Algorithm

Some folks out there find it difficult to understand the Lagrange's algorithm for interpolation polynomial. In this post, I have redeveloped the algorithm in such a format that you will understand and easily implement it in a computer using any of the programming languages of your choice, perhaps java.

The Lagrange's interpolation algorithm is a method or algorithm developed by the French mathematician and astronomer, Joseph Louis Lagrange (1736-1813).

Before we dive into the algorithm that solves the Lagrange form of interpolation polynomial, let us review some concept of polynomial. So that you can trace the footprint that leads to the solution of an interpolation polynomial using pencil and paper. If you are not able to do that, then you will not understand the algorithm. Why? You must know how to solve a problem in other to devise a method (algorithm) that any sane person or device can follow to solve a similar problem.

You will frequently have occasion to estimate intermediate values between precise data points. The most common method used for this purpose is polynomial interpolation. Recall that the general formula for an nth-order polynomial is:

Equation (1)

For n + 1 data points, there is one and only one polynomial of order n that passes through all the points. Polynomial interpolation consists of determining the unique nth-order polynomial that fits n + 1 data points. This polynomial then provides a formula to compute intermediate values.

Although there is one and only one nth-order polynomial that fits n + 1 points, there are a variety of mathematical formats in which this polynomial can be expressed. The two alternatives that are well-suited for computer implementation are: the Newton and the Lagrange polynomials. In this post I will describe the Lagrange polynomials.

Examples of interpolating polynomials: (a) first-order (linear) connecting two points, (b) second-order (quadratic or parabolic) connecting three points, and (c) third-order (cubic) connecting four points. The graph are as show below:


The principle of Lagrange interpolation is that a function f(x) whose values are given at a collection of points is assumed to be approximately represented by a polynomial p(x) that passes through each and every point. The polynomial is called the interpolation polynomial and it is of degree one less than the number of points given. For two data points the interpolating polynomial is taken to be a linear polynomial, as you will see in the below instance. For three data points the interpolating polynomial is taken to be a quadratic, for four data points the interpolating polynomial is taken to be a cubic, and so on.

Equation (2)

Equation (3)

Equation (4)

The Lagrange interpolating polynomial is simply a reformulation of the Newton polynomial that avoids the computation of divided differences. It can also be represented concisely as:

Equation (5)

Equation (6)

Equation (7)

Instance Zero:

From the table below determine/estimate the value of  y, when x = 2.3?

x	1.1	1.7	3.0
F(x) = y	10.6	15.2	20.3

From the table using the Newton method we can say that 2.3 lies between 1.7 and 3.0 in the x field. If this true it therefore means that y lies between 15.2 and 20.3. The Lagrange polynomial equation helps us to determine/solve a problem similar to the one above.

Solution:

From the above equations, we have that;
p(x) = Summation{Li(x)f(xi) = L0(x)f(x0) + L1(x)f(x1) + L2(x)f(x2)
L0(x) = [{(x - x1)(x - x2) / (x0 - x1)(x0 - x2)} * f(x0)]
L0(x) =  [{(2.3 - 1.7)(2.3 - 3.0) / (1.1 - 1.7)(1.1 - 3.0)} * (10.6)]
L0(x) = {(0.6)(-0.7) / (-0.6)(-1.9)} * 10.6 = (-0.42/1.4) * 10.6
L0(x) = -0.36 * 10.6 = -3.90526.
L1(x) = [{(x - x0)(x - x2) / (x1 - x0)(x1 - x2)} * f(x1)]
L1(x) = [{(2.3 - 1.1)(2.3 - 3.0) / (1.7 - 1.1)(1.7 - 3.0)} * (15.2)]
L1(x) = {(1.2)(-0.7) / (0.6)(-1.3)} * 15.2
L1(x) = (-0.84/-0.78) * 15.2
L1(x) = 1.076 * 15.2 = 16.3692.
L2(x) = [{(x - x0)(x - x1) / (x2 - x0)(x2 - x1)} * f(x2)]
L2(x) = [{(2.3 - 1.1)(2.3 - 1.7) / (3.0 - 1.1)(3.0 - 1.7)} * (20.3)}
L2(x) = {(1.2)(0.6) / (1.9)(1.3)} * 20.3
L2(x) = (0.72/2.47) * 20.3
L2(x) = 0.2914 * 20.3 = 5.91740.

From p(x) = Summation{Li(x)f(x) = L0(x)f(x0) + L1(x)f(x1) + L2(x)f(x2)
p(x) = -3.90526 + 16.3692 + 5.91740 = 18.38137.
Therefore When x = 2.3,  f(x) = 18.38137.

Our mission here is not to learn the arithmetic as the one above but to establish a method/steps that solve similar problem at all conditions; at any data points. So that, any sane device or humans can follow the same steps to solve a similar problem at same condition.

We are now worthy to approach the Lagrange's algorithm in its simplest form.

/*
 ************************************************
 * Lagrange Interpolation Polynomial Algorithm  *
 ************************************************
 *
 *************************
 *************************
 * @author Ezehlivinus  **
 * @date Dec 27, 2014   **
 * @time 5:15:16 AM     **
 *************************
 *************************
 */

START
INPUT the number of terms, n
count = 0
count2 = 0
x = 0
y = 0
DO UNTIL count < n
 count = count + 1
 INPUT the Array, arrayx
END DO
DO UNTIL count < n
 count = count + 1
 INPUT the the Arrayy, arrayy
END DO
INPUT the arbitrary value, x

DO UNTIL count < n
 numerator = 1
 denominator = 1
 DO UNTIL count2 < n
  IF count2 != count THEN
   numerator = numerator * (x - arrayx[count2])
   denominator = denominator * (arrayx[count] - arrayx[count2])
  ELSE
   do nothing
  END IF
 END DO
 y = y + (numerator/denominator) * arrayy[count]
END DO

Display x, y
END

Lines 31 to 43 are the lines that do the interpolation computation. Line 31 sets a condition that if true, causes lines 31 to 42 to execute. Otherwise if the condition is false lines 31 to 43 will not execute. numerator and denominator are the variables that will hold the numerator and denominator of the interpolation polynomial equation. Line 34 to 41 causes lines 34 to 41 to execute, if the condition at line 34 is true. If it is false the line 34 to 41 will not execute. Line 35 set up a condition that when true, line 36 and 37 will execute and when false lines 38 to 40 will not execute, before control is transferred to the next line (42), which is under the condition of line 31.

Your question will now be: How does the algorithm or computer know which is x0, x2, x3, ..., and/or xn during execution?

Note that arrayx and arrayy will hold (in group) all the values of x and y respectively. The count and count2 variables determine the position of each element of the arrays. Ones the first element (x0) is entered count = 0, when the second element (x1) is entered count = 1, for third element (x2) count = 2, for fourth element (x3) count = 3, ..., for nth element (xn) count = n. This happens at line 21 - 24. It is similar to y at line 25 to 28.

Instance one:

If x = 2.3, 3, 4, 8 and 7.
Then
x0 = 2.3,
x1 = 3,
x2 = 4,
x3 = 8 and
x4 = 7.

Therefore when:
count = 0, x0 = 2.3;
count = 1, x1 = 3;
count = 2, x2 = 4;
count = 3, x3 = 8;
count = 4, x3 = 7;

Using the above one can easily identify (access or make use of) the elements of x just by calling/specifying their numerical position in the array. So the answer to the above question is that we will number all the elements of x (i.e. the arrayx) call them up using this numbers. The instance below explains further.

Instance two: Using instance one data:
For the first execution (first loop cycle) of the algorithm:
At line 37: arrayx[count] = 2.3,
For the fourth execution of the algorithm:
At line 37: arrayx[count] = 8. How? Because
At the first loop cycle count = 0 = 2.3.
At the second loop cycle count = 1 = 3.
At the third loop cycle count = 2 = 4.
At the Fourth loop cycle count = 3 = 8 and
At the Fifth loop cycle count = 4 = 7.

So apart from the looping characteristics that count variable have it also serve as a place value for the arrays, a number that identifies any element in the array. Note that at line 36 that arrayx[count2] = 1 = 3 = x1. Just to conform with the first bracket (x - x1) in an interpolating polynomial equation . If it is not equal to one (1), but zero (0) as it appears to be during the first execution (loop), the condition at line 35 will not allow line 36 and 37 to execute.

These explanations are also applied to the array arrayy specifically at line 42. All the values of the x and y are stored in their respective arrays.

Since the variables at line 32 and 33 are initialised to be one (1), line 36: "numerator = numerator * (x - arrayx[count2])" is equivalent to the first expression in the numerators of an interpolating polynomial equation as shown below:

Instance three:

1. (x - x1) is equivalent to (x - arrayx[count2]);
2. (x - x1)(x - x2) is equivalent to (x - arrayx[count2])(x - arrayx[count2])
3. (x - x1)(x - x2)(x - x3) is equivalent to (x - arrayx[count2])(x - arrayx[count2])(x - arrayx[count2])
4. (x - x1)(x - x2)(...)(x - xn) equivalent to (x - arrayx[count2])(x - arrayx[count2]) (..)(x - arrayx[count2n]).

How then do we represent any of the expression (numerators only, just for now) in each of the  polynomial after the additional operator '+', to include the second, third and n brackets as indicated in the above expression, in the algorithm. Think before you proceed, I will not be the only one talking, when you have finish sketching your own expression, compare it with the one below. Of course they must not be the same but they should be doing the same thing.

So that - from instance three - we now adjust the four cases as follow:
1. (x - x1) is equivalent to (x - arrayx[count2]);
From equation 2: the numerator of the polynomial: (x - x1) + (x - x0) is equivalent to the numerator: (x - arrayx[count2]) + (x - arrayx[count2]) in the algorithm;

2. (x - x1)(x - x2) is equivalent to (x - arrayx[count2])(x - arrayx[count2]);
From equation 3: The numerator: (x - x1)(x - x2) + (x - x0)(x - x2) + (x - x0)(x - x1), is equivalent to the numerator: (x - arrayx[count2])(x - arrayx[count2]) + (x - arrayx [count2])(x - arrayx[count2]) + (x - arrayx[count2])(x - arrayx[count2]) in the algorithm.

3. (x - x1)(x - x2)(x - x3) is equivalent to (x - arrayx[count2])(x - arrayx[count2])(x - arrayx[count2]);
From equation 4: The numerator: (x - x1)(x - x2)(x - x3) + (x - x0)(x - x2)(x - x3) + (x - x0)(x - x1)(x - x3) + (x - x0)(x - x1)(x - x2), is equivalent to the numerator: (x - arrayx[count2])(x - arrayx[count2])(x - arrayx[count2]) + (x - arrayx[count2])(x - arrayx[count2])(x - arrayx[count2]) + (x - arrayx[count2])(x - arrayx[count2])(x - arrayx[count2]) in the algorithm.

4. (x - x1)(x - x2)(...)(x - xn) equivalent to (x - arrayx[count2])(x - arrayx[count2])(..)(x - arrayx[count2n]);
From equation 5: the numerator: (x - x1)(x - x2)(...)(x - xn) + (x - x0)(x - x2)(...)(x - xn) + (x - x0)(x - x1)(...)(x - xn-1) is equivalent to the numerator: (x - arrayx[count2])(x - arrayx[count2])(...)(x - arrayx[count2n]) + (x - arrayx[count2])(x - arrayx[count2n])(...) + (x - arrayx[count2])(x - arrayx[count2])(...)(x - arrayx[count2n-1]) In the algorithm.

Loop statements helps us to repeat the computation of: (x - arrayx[count2n]) one time in case 1, two times in case 2, three times in case 3 and n times in case 4 for the first, second, third and n brackets in case 1, 2, 3 and 4 respectively. The loop statement that are responsible for this is at line 34 and ends at line 41, if and only if line 34 is true.

When line 34 is false, line 41 stops the execution and transfers control to line 42. It is at line 43 the value of the numerator and the denominator are divided and the quotient is multiplied with arrayy[count]. Where f(x0) = y0 = arrayy[count]. Ones this is done f(x0) in any of the polynomial equation being computed will now have a value which is stored at line 42 in y.

Line 31 specified that count must be less than n before execution will stop. For the sake of this explanation, count is not less n. At the same line, count will be incremented by one. The numerator and denominator variables at line 32 and 33 will be reinitialised to one (1)  for the computation of the second sets of brackets. In a polynomial equation, the second sets of brackets are located after the first sets of bracket followed by the additional operator '+' from left to right. The second set of bracket from case 3 is:

(x - x0)(x - x2)(x - x3)

Operations similar to those performed on the first bracket are performed on the second bracket. Ones the numerator and the denominator are obtain, then the action are transferred to line 42 to obtain y. This is:

y = y + (numerator/denominator) * arrayy[count]

Now to calculate the value of the second y; that is f(x1) = y1 as indicated in the equation. It is equivalent to:

y = y + (numerator/denominator) * arrayy[count], in the algorithm. In a polynomial form the second case (i.e for four data points) it is given as:

p(x1) = (x - x0)(x - x2)(x - x3)/(x1 - x0)(x1 - x2)(x1 - x3) * (fx1).

During execution from the algorithm perspective, this means:
y = (previous value of y) + (new value of the numerator/new value of the denominator) * f(x1).
This process continues for a polynomial passing any given data points.

If you are not happy with  the above algorithm. I bet using your monthly income that you will surely understand this one below:

/*
 ************************************************
 * Lagrange Interpolation Polynomial Algorithm  *
 ************************************************
 *
 *************************
 *************************
 * @author Ezehlivinus  **
 * @date Dec 27, 2014   **
 * @time 5:15:16 AM     **
 *************************
 *************************
 */

START
input the number of terms, n
for (count = 0; count < n; count++)
 Input the array, arrayx
end for
for (count = 0; count < n; count++)
 Input the array, arrayy
end for
input the arbitray value, x
for (count = 0; count < n; count++)
 numerator = 1
 denominator = 1
 for (count2 = 0; count2 < n; count2++)
  if ( count2 != count ) then
   numerator = numerator * (x - arrayx[count2])
   denominator = denominator * (arrayx[count] - arrayx[count2])
  else
   do nothing
  end if
 end for
 y = y + (numerator/denominator) * arrayy[count]
end for
print x, y
END

We are not through yet until you do this one alone:

Question: If line 28 in the second algorithm is changed to:

if ( count2 == count ) then

Then adjust the algorithm where necessary. Make use of the comment box for your answer and/or questions.
At this point you can translate the algorithm into any programming language of your choice. For the java programme that implements the Lagrange’s  interpolation polynomial click here.

SHARE THIS POST WITH YOUR FRIENDS...::

Tweet