Sunday, December 07, 2008

Lagrange's interpretation in python

I've been reading up on Fast Fourier Transform (FFT) from several books on algorithms that I have, TCLR, Tamassia, Sahni, Numerical Recipes, a big set of DSP books. Several of these books focus on interpolation before touching upon the subject of FFT. My favourite book so far for numerical algorithms is R.W. Hamming's Numerical Methods for Scientists and Engineers, Second Edition, Dover publication. Chapter 14, touches upon this topic. I love his explanation and illustration of interpolation versus extrapolation. Thinking along that subject, I decided to do some implementation of my own for interpolation. Sahni (Computer algorithms in C++) has an understandable implementation of the algorithm. The biggest drawback was implementing my own polynomial class. This is where, I think Python scores. I found that SciPy implements a polynomial class.

Here is the code for interpolation

def lagrange(x):
tmp = scipy.poly1d([0])

for i in x.keys():
denom = 1.0
for j in x.keys():
if (i != j):
tmp = scipy.poly1d([1,-j])
numerator = numerator * tmp
denom = denom * (i - j)
tmp = (numerator/denom) * x.get(i)
result = result + tmp

return result

input = {0:5, 1:10, 2:21}
print lagrange(input)
NOTE: Blogger keeps messing up my indentation, I am yet to figure out why.

Sahni's input for the example is {0:1, 1:10, 2:21} and I was surprised to see my program return a different output. Fixing the input, showed the correct results. That reminds me, somebody please ask Sahni to start maintaining useful errata for his books. I searched on the web and found nothing useful.

Running the same algorithm on y = log(x) and input = {1:0, 2:0.3010, 3:0.4771, 4:0.6021}, gave me 00.0123 x^3 - 0.1362 x^2 + 0.6236 x - 0.4997, which corresponds to the cubical approximation of log x (R W Hamming, page 233).


Anonymous said...
This comment has been removed by a blog administrator.
Mollier said...

Hi, I just wrote a blog-post about Lagrange polynomials and I hope you don't mind that I modified your code for my use.
Thanks for sharing.

Mollier said...

Sorry, I did not post a link to my blog so that you can approve.

Balbir Singh said...

Hey, Mollier

No issues, please feel free to use it. I looked at your blog entry it is quite cool!

Dynamic programming for the binomial coefficient

More fun things, this time with some visualisation of what happens when memoisation is used and what happens when we don't. I don'...