The Lagrange interpolation formula can be programmed in Mathematica. The idea is to turn a list of elements into a product of differences excluding a specified dropped element. Here is how it can be done using some Mathematica shorthand.
data = {a, b, c, d, e}
Drop[data, {4}]
x - %
Times @@ %
We can create a single function to do this, where
"data" is the list (array) of elements,
"x" is an element the other elements are to be subtracted from, and
"drop" is the integer index of the location of the element in "data" to be dropped.
L[data_,x_,drop_] := Times @@ (x - Drop[data, {drop}])
Let's try it out.
![[Graphics:Images/index_gr_1.gif]](Images/index_gr_1.gif)
Below I have written the base code for a function lip[] that returns a Lagrange Interpolating Polynomial, where
"a" is a list of abscissas,
, (i = 0, 1, ..., n),
"b" is a list of ordinates,
, (i = 0, 1, ..., n), and
"p" is the name of the independent variable
If the abscissas are not distinct, the function will not do any computation, and will return an error message.
Here is a sample call to lip[] :
lip[{0, 1, 4}, {1, 2, 17}, x]
This particular call would return the Lagrange interpolating polynomial in x for the points (0, 1), (1, 2), (4, 17).
I have left part of the code blank. We will fill it in during lab.
lip::nds = "The abscissas `1` are not distinct";
lip[a_,b_,p_] := If [Sort[a] != Union[a],
Message[lip::nds,a],
len=Length[a];
poly = 0;
For[i = 1, i <= len, i++, poly = poly + (what goes here?) ];
poly
]
Let's try out our new function.
lip[{0, 1, 4}, {1, 2, 17}, x]
The function InterpolatingPolynomial[] in Mathematica gives the Newton form of the interpolating polynomial to a set of data.
This call returns the interpolating polynomial in the variable x that passes through the points {x0,f0}, {x1,f1}, {x2,f2} :
![[Graphics:Images/index_gr_4.gif]](Images/index_gr_4.gif)
This call defines the interpolating polynomial ip in the variable x that passes through the points {x0,f0}, {x0+h,f1}, {x0 + 2h,f2} :
![[Graphics:Images/index_gr_5.gif]](Images/index_gr_5.gif)
Evaluating ip at x = x0 + p*h :
![[Graphics:Images/index_gr_6.gif]](Images/index_gr_6.gif)
Let's simplify the result. This is the Newton Forward Difference Formula!
![[Graphics:Images/index_gr_7.gif]](Images/index_gr_7.gif)
Let's try the InterpolatingPolynomial[] function out on our previous data set.
![[Graphics:Images/index_gr_8.gif]](Images/index_gr_8.gif)
That doesn't look like the result we got from determining the Lagrange form! What's wrong?!?!?!
Let's simplify the result from the Lagrange form and check again.
![[Graphics:Images/index_gr_9.gif]](Images/index_gr_9.gif)
Here we present another method of interpolating a set of points. The method uses undetermined coefficients. The idea is to fit the set of n+1 points with a polynomial of degree n where the polynomial has the form:
+
+ ... +
+
= 0
This method gives a list of n+1 equations with n+1 unknowns, which can be solved by Gaussian elimination.
Function muc takes a list of abscissas,
, (i = 0, 1, ..., n), a list of ordinates,
, (i = 0, 1, ..., n), and a variable p, and returns the interpolating polynomial calculated by the method of undetermined coefficients (determining the coefficients of the interpolating polynomial by solving a set of (n + 1) linear equations in (n + 1) unknowns - the unknowns are the coefficients). If the abscissas are not distinct, the function will not do any computation, and will return an error message.
muc::nds = "The abscissas `1` are not distinct";
muc[a_,b_,p_] := If [Sort[a] != Union[a],
Message[muc::nds,a],
len=Length[a];
LinearSolve[ Table[
Join[{1}, a[[i]]^Range[1, len - 1]], {i, len}], b]
. (p^Range[0, len - 1])
]
The interpolating polynomial for a given set of data is unique. This can be illustrated nicely by generating interpolating polynomials using three different methods, and then plotting them all together.
Lform = lip[{-1, 1, 3, 7}, {4, 1, 2, 9}, x]
Nform = InterpolatingPolynomial[{{-1,4},{1,1},{3,2},{7,9}},x]
MUCform = muc[{-1, 1, 3, 7}, {4, 1, 2, 9}, x]
![[Graphics:Images/index_gr_16.gif]](Images/index_gr_16.gif)