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.

Start with defining

\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.

About these ads

One comment

  1. Natalie Schneider

    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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 44 other followers

%d bloggers like this: