Lagrange interpolation in python

lagrange.py

import numpy as np
import matplotlib.pyplot as plt
import sys

def main():
	if len(sys.argv) == 1 or "-h" in sys.argv or "--help" in sys.argv:
		print "python lagrange.py <x1.y1> .. <x_k.y_k>"
		print "Example:"
		print "python lagrange.py 0.1 2.4 4.5 3.2"
		exit()
	points = []
	for i in xrange(len(sys.argv)): 
		if i != 0: 
			points.append((int(sys.argv[i].split(".")[0]),int(sys.argv[i].split(".")[1]))) 

	#points =[(0,0),(25,30),(50,10), (57,0)]
	P = lagrange(points)
	nr = 2
	print "(" + str(points[nr][0]) + ", " + str(points[nr][1]) +") P(" + str(points[nr][0]) +")= " +str(P(points[nr][0]))
	plot(P, points)

def plot(f, points):
	x = range(-10, 100)
	y = map(f, x)
	print y
	plt.plot( x, y, linewidth=2.0)
	x_list = []
	y_list = []
	for x_p, y_p in points:
		x_list.append(x_p)
		y_list.append(y_p)
	print x_list
	print y_list
	plt.plot(x_list, y_list,  'ro')

	plt.show()


def lagrange(points):
	def P(x):
		total = 0
		n = len(points)
		for i in xrange(n):
			xi, yi = points[i]
			
			def g(i, n):
				
				tot_mul = 1
				for j in xrange(n):
					if i == j:
						continue
					xj, yj = points[j]
					tot_mul *= (x - xj) / float(xi - xj)
					
				return tot_mul 

			total += yi * g(i, n)
		return total
	return P 


if __name__ == "__main__":
	main()

Hello,

Would it be possible to get the literal expression of polynom P?

I mean, printing out something like "P(x) = 5x^3-2x^2+6x+7".

It would be another method with x being a symbolic variable. See my version : https://gist.github.com/aurelienpierre/1d9826e7db078e048bf437e516a7a4b2

Your 1D lagrange is neater than the version published by Stoecker. He gives a rough suggestion at a 2D lagrange but I'm not winning with the nested loops. Have you tried 2D?

Here's my take on a 2D version.

"""
Lagrange 2D
Stoecker "Design of Thermal Systems", 2nd ed. page 62.
Bruce Wernick
22 September 2017 5:21:48
"""

from __future__ import division

def lagrange2(x, y, X, Y, Z):
  m = len(X)
  n = len(Y)
  z = 0.0
  for i in range(m):
    for j in range(n):
      p, q = 1.0, 1.0
      for k in range(m):
        if i == k: continue
        p *= (x - X[k])
        q *= (X[i] - X[k])
      for k in range(n):
        if j == k: continue
        p *= (y - Y[k])
        q *= (Y[j] - Y[k])
      z += Z[j][i] * p / q
  return z

X = [20,25,30,35]
Y = [5,10,15]
Z = [[100,120,140,160], [130,150,170,190], [160,180,200,220]]

print lagrange2(40, 10, X, Y, Z)

It seems to mess up the spacing in the code?

Is this code restricted by the order of the polynomial that it outputs? ie, will it was always output n degree polynomials, or is it variable, and will depend on the input?