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])
result=scipy.poly1d([0])

for i in x.keys():
numerator=scipy.poly1d([1])
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).
Post a Comment