# coding: utf-8 # # Orthogonal polynomials # In[33]: import numpy as np import numpy.linalg as la import matplotlib.pyplot as pt # ## Mini-Introduction to `sympy` # In[34]: import sympy as sym # Enable "pretty-printing" in IPython sym.init_printing() # Make a new `Symbol` and work with it: # In[35]: x = sym.Symbol("x") myexpr = (x**2-3)**2 myexpr # In[36]: myexpr = (x**2-3)**2 myexpr myexpr.expand() # In[37]: sym.integrate(myexpr, x) # In[38]: sym.integrate(myexpr, (x, -1, 1)) # ## Orthogonal polynomials # Now write a function `inner_product(f, g)`: # In[39]: def inner_product(f, g): return sym.integrate(f*g, (x, -1, 1)) # Show that it works: # In[40]: inner_product(1, 1) # In[41]: inner_product(1, x) # Next, define a `basis` consisting of a few monomials: # In[42]: #basis = [1, x, x**2, x**3] basis = [1, x, x**2, x**3, x**4, x**5] # And run Gram-Schmidt on it: # In[43]: orth_basis = [] for q in basis: for prev_q in orth_basis: q = q - inner_product(prev_q, q)*prev_q q = q / sym.sqrt(inner_product(q, q)) orth_basis.append(q) # In[44]: orth_basis # These are called the *Legendre polynomials*. # -------------------- # What do they look like? # In[45]: mesh = np.linspace(-1, 1, 100) pt.figure(figsize=(8,8)) for f in orth_basis: f = sym.lambdify(x, f) pt.plot(mesh, [f(xi) for xi in mesh]) # ----- # These functions are important enough to be included in `scipy.special` as `eval_legendre`: # # **!!** Careful: The Scipy versions do not have norm 1. # # In[62]: import scipy.special as sps for i in range(10): pt.plot(mesh, sps.eval_legendre(i, mesh)) # What can we find out about the conditioning of the generalized Vandermonde matrix for Legendre polynomials? # In[69]: n = 20 xs = np.linspace(-1, 1, n) V = np.array([ sps.eval_legendre(i, xs) for i in range(n) ]).T la.cond(V)