# Orthogonal polynomials


```python
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt
```

## Mini-Introduction to `sympy`


```python
import sympy as sym

# Enable "pretty-printing" in IPython
sym.init_printing()
```

Make a new `Symbol` and work with it:


```python
x = sym.Symbol("x")

myexpr = (x**2-3)**2
myexpr
```


```python
myexpr = (x**2-3)**2
myexpr
myexpr.expand()
```


```python
sym.integrate(myexpr, x)
```


```python
sym.integrate(myexpr, (x, -1, 1))
```

## Orthogonal polynomials

Now write a function `inner_product(f, g)`:


```python
def inner_product(f, g):
    return sym.integrate(f*g, (x, -1, 1))
```

Show that it works:


```python
inner_product(1, 1)
```


```python
inner_product(1, x)
```

Next, define a `basis` consisting of a few monomials:


```python
basis = [1, x, x**2, x**3]
#basis = [1, x, x**2, x**3, x**4, x**5]
```

And run Gram-Schmidt on it:


```python
orth_basis = []

for q in basis:
    for prev_q in orth_basis:
        q = q - inner_product(prev_q, q)*prev_q / inner_product(prev_q,prev_q)
    orth_basis.append(q)

legendre_basis = [orth_basis[0],]

#to compute Legendre polynomials need to normalize so that q(1)=1 rather than ||q||=1
for q in orth_basis[1:]:
    q = q / q.subs(x,1)
    legendre_basis.append(q)
```


```python
legendre_basis
```

These are called the *Legendre polynomials*.

--------------------
What do they look like?


```python
mesh = np.linspace(-1, 1, 100)

pt.figure(figsize=(8,8))
for f in legendre_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`:


```python
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?


```python
#keep
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)
```

The Chebyshev basis can similarly be defined by Gram-Schmidt, but now with respect to a different inner-product weight function,
$$w(x) = 1/\sqrt{1-x^2}.$$


```python
w = 1 / sym.sqrt(1-x**2)
def cheb_inner_product(f, g):
    return sym.integrate(w*f*g, (x, -1, 1))

orth_basis = []

for q in basis:
    for prev_q in orth_basis:
        q = q - cheb_inner_product(prev_q, q)*prev_q / cheb_inner_product(prev_q,prev_q)
    orth_basis.append(q)

cheb_basis = [1,]

#to compute Legendre polynomials need to normalize so that q(1)=1 rather than ||q||=1
for q in orth_basis[1:]:
    q = q / q.subs(x,1)
    cheb_basis.append(q)
cheb_basis
```


```python
for i in range(10):
    pt.plot(mesh, np.cos(i*np.arccos(mesh)))
```

Chebyshev polynomials achieve similar good, but imperfect conditioning on a uniform grid (but perfect conditioning on a grid of Chebyshev nodes).


```python
#keep
n = 20
xs = np.linspace(-1, 1, n)
V = np.array([
    np.cos(i*np.arccos(xs))
    for i in range(n)
]).T

la.cond(V)
```


```python

```

14d10b36ecdbd06408cc2c8cfc59f098e5c9b2be
ipynb
Jupyter Notebook
interpolation/Orthogonal Polynomials.ipynb
JiaheXu/MATH
9cb2b412ba019794702cacf213471742745d17a6
MIT
interpolation/Orthogonal Polynomials.ipynb
JiaheXu/MATH
9cb2b412ba019794702cacf213471742745d17a6
MIT
interpolation/Orthogonal Polynomials.ipynb
JiaheXu/MATH
9cb2b412ba019794702cacf213471742745d17a6
MIT
Qwen/Qwen-72B
1. YES
2. YES