# Import symbolics and plotting functions.
from pylab import plot, savefig, grid, legend, clf
from numpy import arange
from sympy import (var, pi, factorial, sin, cos, log, exp, tan, atan, asin,
acos, sinh, cosh, tanh, asinh, acosh, lambdify)
# Use the symbol 'x' for independent variable.
var("x")
# Basic function evaluating at point 'x'
# the Lagrange polynomial given via a list of
# [x, y] pairs.
def lagrange(x, list):
r = 0
for pi in list:
val = 1
for pj in list:
if pi != pj:
val *= (float(x) - pj[0])/(pi[0]-pj[0])
r += val*pi[1]
return r
# Case 1: Interpolation points are defined via
# arbitrary arrays of x- and y-coordinates.
def lagrange_1(a, b, x_list, y_list, n=100):
list = zip(x_list, y_list)
h = (b-a)/float(n)
x_array = arange(a, b+h, h)
lagr_poly = [lagrange(xx, list) for xx in x_array]
clf()
label = "Lagrange polynomial L(x)"
plot(x_array, lagr_poly, "b-", label=label)
plot(x_list, y_list, "ko")
legend()
grid(True)
savefig("a.png")
# Case 2: x-coordinates of interpolation points are
# defined via an arbitrary array x_list, y-coordinates
# via function values f(x).
def lagrange_2(a, b, x_list, f, x, n=100):
f = lambdify(x, f, modules=["math"])
y_list = [f(xx) for xx in x_list]
list = zip(x_list, y_list)
h = (b-a)/float(n)
x_array = arange(a, b+h, h)
lagr_poly = [lagrange(xx, list) for xx in x_array]
f_graph = [f(xx) for xx in x_array]
clf()
label = "Function f(x)"
plot(x_array, f_graph, "g-", label=label)
plot(x_list, y_list, "ko")
label = "Lagrange polynomial L(x)"
plot(x_array, lagr_poly, "b-", label=label)
legend()
grid(True)
savefig("a.png")
clf()
error_graph = [float(f(xx) - lagrange(xx, list)) for xx in x_array]
label = "Error E(x) = f(x) - L(x)"
plot(x_array, error_graph, "b-", label=label)
legend()
grid(True)
savefig("b.png")
# calculate approximate error maximum (over plotting points)
err_abs = [abs(float(f(xx) - lagrange(xx, list))) for xx in x_array]
print "Max error:", max(err_abs)
# Case 3: interpolation points are equidistributed
# in the interval [a,b], y-coordinates are defined
# via a function f(x).
def lagrange_3(a, b, num_pts, f, x, n=100):
f = lambdify(x, f, modules=["math"])
num_elem = num_pts-1
h0 = (b-a)/float(num_elem)
x_list = [float(a + i*h0) for i in range(0, num_pts)]
# demonstration of inoptimality of equidistant points:
# move points out of center and see error decrease
#x_list[1] -= 0.1
#x_list[-2] += 0.1
y_list = [f(xx) for xx in x_list]
list = zip(x_list, y_list)
h = (b-a)/float(n)
x_array = arange(a, b+h, h)
lagr_poly = [lagrange(xx, list) for xx in x_array]
f_graph = [f(xx) for xx in x_array]
clf()
label = "Function f(x)"
plot(x_array, f_graph, "g-", label=label)
label = "Lagrange polynomial L(x)"
plot(x_array, lagr_poly, "b-", label=label)
plot(x_list, y_list, "ko")
legend()
grid(True)
savefig("a.png")
clf()
error_graph = [f(xx) - lagrange(xx, list) for xx in x_array]
label = "Error E(x) = f(x) - L(x)"
plot(x_array, error_graph, "b-", label=label)
legend()
grid(True)
savefig("b.png")
# calculate approximate error maximum (over plotting points)
err_abs = [abs(float(f(xx) - lagrange(xx, list))) for xx in x_array]
print "Max error:", max(err_abs)
# Case 4: x-coordinates of interpolation points are
# the Chebyshev points, y-coordinates are defined
# via function f(x).
def lagrange_4(a, b, num_pts, f, x, n=100):
f = lambdify(x, f, modules=["math"])
num_elem = num_pts - 1
x_list = [float(cos(i*pi/num_elem)*(b-a)/2) for i in range(0, num_pts)]
# demonstration of optimality of Chebyshev points:
# move any point and see error increase
#x_list[2] += 0.2
#x_list[-3] -= 0.2
y_list = [f(xx) for xx in x_list]
list = zip(x_list, y_list)
h = (b-a)/float(n)
x_array = arange(a, b+h, h)
lagr_poly = [lagrange(xx, list) for xx in x_array]
f_graph = [f(xx) for xx in x_array]
clf()
label = "Function f(x)"
plot(x_array, f_graph, "g-", label=label)
label = "Lagrange polynomial L(x)"
plot(x_array, lagr_poly, "b-", label=label)
plot(x_list, y_list, "ko")
legend()
grid(True)
savefig("a.png")
clf()
error_graph = [f(xx) - lagrange(xx, list) for xx in x_array]
label = "Error E(x) = f(x) - L(x)"
plot(x_array, error_graph, "b-", label=label)
legend()
grid(True)
savefig("b.png")
# calculate approximate error maximum (over plotting points)
err_abs = [abs(float(f(xx) - lagrange(xx, list))) for xx in x_array]
print "Max error:", max(err_abs)
