I want to interpolate a polynomial with the Lagrange method, but this code doesn't work:
def interpolate(x_values, y_values):
def _basis(j):
p = [(x - x_values[m])/(x_values[j] - x_values[m]) for m in xrange(k + 1) if m != j]
return reduce(operator.mul, p)
assert len(x_values) != 0 and (len(x_values) == len(y_values)), 'x and y cannot be empty and must have the same length'
k = len(x_values)
return sum(_basis(j) for j in xrange(k))
I followed Wikipedia, but when I run it I receive an IndexError at line 3!
Thanks
Try
def interpolate(x, x_values, y_values):
def _basis(j):
p = [(x - x_values[m])/(x_values[j] - x_values[m]) for m in xrange(k) if m != j]
return reduce(operator.mul, p)
assert len(x_values) != 0 and (len(x_values) == len(y_values)), 'x and y cannot be empty and must have the same length'
k = len(x_values)
return sum(_basis(j)*y_values[j] for j in xrange(k))
You can confirm it as follows:
>>> interpolate(1,[1,2,4],[1,0,2])
1.0
>>> interpolate(2,[1,2,4],[1,0,2])
0.0
>>> interpolate(4,[1,2,4],[1,0,2])
2.0
>>> interpolate(3,[1,2,4],[1,0,2])
0.33333333333333331
So the result is the interpolated value based on the polynomial that goes through the points given. In this case, the 3 points define a parabola and the first 3 tests show that the stated y_value is returned for the given x_value.