## Easy Lagrangian interpolation

As everybody knows, the formula for Lagrangian interpolation though points $(x_i,y_i)$, where all the $x_i$ values are distinct, is given by

$\displaystyle{P(x)=\sum_{i=0}^n\left(\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\right)y_i}$

It is easy to see that this works, as the product is equal to zero if $x=x_k\ne x_i$ and is equal to one if $x=x_i$.

This isn’t particularly easy to program, and in fact methods such as Newton’s divided differences or Neville’s algorithm are used to obtain the polynomial.

But a method introduced by Gabor Szegö provides a very simple form of Lagrange’s polynomial.

$\pi(x)=(x-x_0)(x-x_1)\cdots(x-x_n)$.

Then

$\displaystyle{\frac{d}{dx}\pi(x)=\sum_{i=0}^n\Biggl(\prod_{\substack{0\le j\le n\\j\ne i}}(x-x_j)\Biggr)}.$

This is obtained immediately from the product rule for arbitrary products, and using the fact that the derivative of each individual term of $\pi(x)$ is one. This means that in particular

$\displaystyle{\pi'(x_i) = \prod_{\substack{0\le j\le n\\j\ne i}}(x_i-x_j)}$

because all other terms of $\pi'(x)$ contain an $(x-x_i)$ term and therefore produce zero for $x=x_i$.

This means that the coefficient of $y_i$ in the Lagrangian polynomial can be written as

$\displaystyle{\frac{\pi(x)}{(x-x_i)\pi'(x_i)}}$

and so the polynomial can be written as

$\displaystyle{P(x)=\sum_{i=0}^n\frac{\pi(x)}{(x-x_i)\pi'(x_i)}y_i}$.

Here is an example, with $x=-2,1,3,4$ and $y=-3,-27,-23,3$. First, Maxima:

(%o1) xs:[-2,1,3,4]; (%o2) ys:[-3,-27,-23,3]; (%o3) define(pi(x),product(x-xs[i],i,1,4)); (%o4) define(pid(x),diff(pi(x),x)); (%o5) define(P(x),sum(pi(x)/(x-xs[i])/pid(xs[i])*ys[i],i,1,4)); (%o6) ratsimp(P(x));

$x^3-11x-17$

Now Sage:

sage: xs = [-2,1,3,4]
sage: ys = [-3, -27, -23, 3]
sage: pi(x) = prod(x-i for i in xs); pi(x)
(x - 4)*(x - 3)*(x - 1)*(x + 2)
sage: pid(x) = diff(pi(x),x);
sage: P(x) = sum(pi(x)/(x-i)/pid(i)*j for (i,j) in zip(xs,ys))
sage: P(x).collect(x)
x^3 - 11*x - 17


And finally, on the TI-nspire CAS calculator:

$n:=4$

$xs:=\{-2,1,3,4\}$

$ys:=\{-3,-27,-23,3\}$

$\displaystyle{{\rm Define}\quad q(x)=\prod_{i=1}^n(x-xs[i])}$

$\displaystyle{qs:={\rm seq}\left(\left.\frac{d}{dx}(q(x))\right|x=xs[i],i,1,n\right)}$

$\displaystyle{{\rm Define }\quad p(x)=\sum_{i=1}^n\left(\frac{q(x)}{(x-xs[i]) \cdot qs[i]}ys[i]\right)}$

$p(x)\hspace{20mm} x^3-11\cdot x-17$

Note that all the mathematical symbols and programming constructs are available through the calculator’s menus.

Notice that this approach requires no programming with loops or conditions at all – it’s all done with simple sums and products.

### One Response

1. Placing this in a program editor would make it even easier to use in the Nspire CAS. Thanks for the syntax of how to do it though.

use instead a program editor, with the two inputs as the lists

also use n:= count(xs) to determine the length rather than typing a new length every time.

for some reason when I used the input variables directly it wouldn’t work, so I redefined them and it did. Maybe you can make it better, but here is the code for that program which will work on any two lists of the same length (Crashes if y is shorter)

-code

Define laginterp(xls,yls)=
Prgm
:xs:=xls
:ys:=yls
:n:=count(xs)
:Define q(x)=∏(x-xs[i],i,1,n)
:qs:=seq(d/dx(q(x),x)|x=xs[i],i,1,n)
:Define p(x)=∑(((q(x))/((x-xs[i])*qs[i]))*ys[i],i,1,n)
:Disp p(x)
:EndPrgm

The AIM Network

The Australian Independent Media Network