# Import symbolics, plotting, linear algebra, Legendre
# polynomials, and numerical quadrature functions.
from pylab import plot, savefig, grid, legend, clf, pcolor, spy, axis
from numpy import arange, zeros
from math import (pi, factorial, sin, cos, log, exp, tan, atan, asin,
acos, sinh, cosh, tanh, asinh, acosh, sqrt)
from scipy.linalg import solve
from scipy.special.orthogonal import legendre
from scipy.special.orthogonal import p_roots
from scipy.integrate import quadrature
# Define Legendre basis.
def P(i, x):
return legendre(i)(x)
# The Legendre polynomials are defined in the interval (-1,1)
# while the function f(x) in the interval (a, b). Therefore we
# need a reference mapping from (-1, 1) onto (a, b)
def refmap(a, b, x):
return (a+b)/2. + x*(b-a)/2.
# Evaluate at 'x' linear combination of Legendre polynomials
# with the coefficient vector 'vec'
def lin_comb(x, vec):
val = 0
n = len(vec)
for i in range(0, n):
val = val + vec[i]*P(i, x)
return val
# Integrate product of f(x) with fn(i)
def gauss_quad(a, b, f, fn, i, quad_x, quad_w):
val = 0
for j in range(0, len(quad_x)):
val += f(refmap(a, b, quad_x[j])) * fn(i, quad_x[j]) * quad_w[j]
norm_squared = 2./(2.*i + 1.)
return val/norm_squared
# Calculate L2 error norm
def calculate_L2_error(a, b, f, rhs, quad_x, quad_w):
err = 0
for j in range(0, len(quad_x)):
val = f(refmap(a, b, quad_x[j])) - lin_comb(quad_x[j], rhs)
err += val**2 * quad_w[j]
return sqrt(err)
# Calculate projection coefficients. Use the fact that
# after normalization, Legendre polynomials are orthogonal
# on (-1,1) in the L2-product (and thus no matrix inversion
# is necessary).
def calc_projection_coeffs(a, b, deg, f, num_gauss_points, print_rhs):
n = deg + 1
rhs = zeros(n) # zero vector
# preparing Gauss quadrature
print "Preparing Gauss quadrature..."
quad_x, quad_w = p_roots(num_gauss_points)
print "Done."
print "num_gauss_points =", num_gauss_points
print "WARNING: You should experiment with num_gauss_points"
print " and increase it until results converge."
# fill the rhs vector
print "Calculating integrals..."
for i in range(0, deg+1): # loop over Legendre polynomials
# calculating the product of P(i) with f
rhs[i] = gauss_quad(a, b, f, P, i, quad_x, quad_w)
print "Done."
# printing right-hand side rhs
if print_rhs: print "Projection coefficients:", rhs
# calculate L2 error numerically
print "Calculating L2 error..."
err_L2 = calculate_L2_error(a, b, f, rhs, quad_x, quad_w)
print "L2 error:", err_L2
return rhs
# Plot the projection result, which is a linear combination of
# Legendre polynomials with the coefficients from the vector vec.
def plot_result(a, b, vec, f, m=200):
print "Plotting result..."
h_ref = 2./m # 2 is the length of the reference interval (-1, 1)
ref_array = arange(-1, 1+h_ref, h_ref)
x_array = [refmap(a, b, xx) for xx in ref_array]
y_array = [lin_comb(xx, vec) for xx in ref_array]
f_array = [f(refmap(a, b, xx)) for xx in ref_array]
e_array = [f(refmap(a, b, xx)) - lin_comb(xx, vec) for xx in ref_array]
clf()
label = "Approximation P(x)"
axis('equal')
plot(x_array, y_array, "b-", label=label)
label = "Function f(x)"
plot(x_array, f_array, "g-", label=label)
legend()
grid(True)
savefig("a.png")
clf()
label = "Error E(x) = f(x) - P(x)"
plot(x_array, e_array, "b-", label=label)
legend()
grid(True)
savefig("err.png")
print "Done."
# approximate error maximum (over plotting points)
print "Calculating max error..."
err_abs = [abs(f(refmap(a, b, xx)) - lin_comb(xx, vec)) for xx in ref_array]
err_max = max(err_abs)
print "Max error:", err_max
return
# Plot the basis.
def plot_basis(a, b, deg, m=200):
print "Plotting basis functions..."
h_ref = 2./m # 2 is the length of the reference interval (-1, 1)
x_array = arange(-1, 1+h_ref, h_ref)
for i in range(0, deg+1):
clf()
label = "Legendre polynomial P%d(x)" % (i)
y_array = [P(i, xx) for xx in x_array]
plot(x_array, y_array, "b-", label=label)
legend()
grid(True)
savefig("leg%d.png" % (i))
print "Done."
return
# Construct and plot the best approximation P(x) of a given function f(x)
# in the space of polynomials of degree less than or equal to 'deg',
# defined in interval (a,b).
def least_squares(a, b, f, deg, plot_basis_fns=False, print_rhs=True, m=200):
a = float(a)
b = float(b)
deg = int(deg)
# WARNING: The number of Gauss points below is chosen heuristically.
# You should experiment with it and increase it if needed.
num_gauss_points = deg + 10 # remember: quadrature order = 2*this - 1
# Perform orthogonal L2-projection of f(x) onto the linear
# space generated by the Legendre polynomials up to degree deg.
vec = calc_projection_coeffs(a, b, deg, f, num_gauss_points, print_rhs)
# plot the approximation P(x) along with f(x) in interval (a, b)
plot_result(a, b, vec, f)
# Plot Legendre polynomials up to degree deg in interval (-1,1).
if plot_basis_fns: plot_basis(a, b, deg)
return
