# M2AA3 Chapter 2, Lesson 2 - Generating Orthogonal Polynomials


```python
#Setup for Sympy
import sympy as sp
```

We wish to construct a list of polynomial $\phi_n(x)$ which is orthogonal with respect to the inner product
$$\langle f(x),g(x) \rangle = \int_a^b w(x)f(x)g(x)dx$$
where $w(x) > 0$ is the weight function. Let $\|f(x)\|^2 = \langle f(x),f(x) \rangle$. We may use the following recurrence relation:
$$\phi_j(x) = \bigg(x - \frac{\langle x\phi_{j-1}(x),\phi_{j-1}(x) \rangle}{\|\phi_{j-1}(x)\|^2} \bigg)\phi_{j-1}(x) - \frac{\|\phi_{j-1}(x)\|^2}{\|\phi_{j-2}(x)\|^2}\phi_{j-2}(x), j \geq 1$$ 
where $\phi_{-1}(x)=0$ and $\phi_{0}(x)=1$. Therefore we have
$$\phi_{1}(x) = x - \frac{\langle x,1 \rangle}{\| 1 \|^2} = x - \frac{\int_a^b x w(x)dx}{\int_a^b w(x)dx}$$
Our aim is to define a function 'ortho', which inputs $w(x)$, $a$, $b$ and $r$ and output a list $\phi_0(x), ..., \phi_r(x)$.

Step 1 - Define inner product and norm


```python
x = sp.symbols('x')
w = sp.Function('w')
f = sp.Function('f')
g = sp.Function('g')
```


```python
def inner(w,a,b,f,g):
    output = sp.integrate(w*f*g, (x,a,b))
    return output

def norm(w,a,b,f):
    return sp.sqrt(inner(w,a,b,f,f))
```

Step 2 - Write the recursive function


```python
phi = sp.Function('phi')
phiminus1 = sp.Function('phiminus1')
phiminus2 = sp.Function('phiminus2')
```


```python
def ortho(w,a,b,r):
    phiminus2 = 1
    phiminus1 = x - inner(w,a,b,x,1)/(norm(w,a,b,1))**2
    if r == 0:
        return [phiminus2]
    elif r == 1:
        return [phiminus2,phiminus1]
    else:
        philist = [phiminus2,phiminus1]
        for i in range(r-1):
            phi = (x - inner(w,a,b,x*phiminus1,phiminus1)/(norm(w,a,b,phiminus1))**2)*phiminus1 - ((norm(w,a,b,phiminus1)/norm(w,a,b,phiminus2))**2)*phiminus2
            phi = sp.simplify(phi)
            philist.append(phi)
            phiminus2 = phiminus1
            phiminus1 = phi
    return philist
```

## Legendre Polynomials
Here is an example of list of polynomials which are orthogonal with $w(x) = 1$. These are called Legendre Polynomials.


```python
olist = ortho(1,-1,1,5)
print(olist)
```

    [1, x, x**2 - 1/3, x*(x**2 - 3/5), x**4 - 6*x**2/7 + 3/35, x*(63*x**4 - 70*x**2 + 15)/63]


Observe that they are very similar to the Chebyshev Polynomials. Recall that they are defined as followed:
$$L_j(x) = \frac{2j-1}{j} x L_{j-1}(x) - \frac{j-1}{j} L_{j-2}(x)$$
where $L_0(x) = 1, L_1(x) = x$. If $\phi_j(x)$ is the orthogonal polynomials obtained, then 
$\phi_j(x) = \frac{L_j(x)}{a_j}, j \geq 1$, where $a_j$ is the leading coefficient of $L_j(x)$.

To verify this let's generate the first $r$ Legendre Polynomials.


```python
P = sp.Function('P')
Pminus1 = sp.Function('Pminus1')
Pminus2 = sp.Function('Pminus2')
```


```python
def Leb(r):
    Pminus2 = 1
    Pminus1 = x
    if r == 0:
        return [Pminus2]
    elif r == 1:
        return [Pminus2,Pminus1]
    else:
        Plist = [Pminus2,Pminus1]
        for i in range(2,r+1):
            P = (sp.Rational(2*i-1,i))*x*Pminus1 - sp.Rational(i-1,i)*Pminus2
            P = sp.simplify(P)
            Plist.append(P)
            Pminus2 = Pminus1
            Pminus1 = P
    return Plist
```

Then we verify that the polynomials are the same.


```python
Leblist = Leb(3)
print(Leblist)
Lebdivided = [1]+[sp.simplify(Leblist[i]/sp.polys.polytools.LC(Leblist[i])) for i in range(1,len(Leblist))]
print(Lebdivided)
plist = ortho(1,-1,1,3)
print(plist)
print(plist == Lebdivided)
```

    [1, x, 3*x**2/2 - 1/2, x*(5*x**2 - 3)/2]
    [1, x, x**2 - 1/3, x*(x**2 - 3/5)]
    [1, x, x**2 - 1/3, x*(x**2 - 3/5)]
    True


We may find the roots of polynomials.


```python
sp.solve(olist[3],x)
```




    [0, -sqrt(15)/5, sqrt(15)/5]



## Chebyshev Polynomials
Here is another example of list of polynomials which are orthogonal with $w(x) = \frac{1}{\sqrt{1 - x^2}}$.


```python
tlist = ortho(1/sp.sqrt(1-x**2),-1,1,3)
print(tlist)
```

    [1, x, x**2 - 1/2, x*(x**2 - 3/4)]


Observe that they are very similar to the Chebyshev Polynomials. Recall that they are defined as followed:
$$T_j(x) = 2x T_{j-1}(x) - T_{j-2}(x), j \geq 2$$
where $T_0(x) = 1, T_1(x) = x$. It could be proved that $T_j = \cos (j \cos^{-1} x)$.
If $\phi_j(x)$ is the orthogonal polynomials obtained, then $\phi_j(x) = \frac{T_j(x)}{2^{j-1}}, j \geq 1$. <br>

To verify this let's generate the first $r$ Chebyshev Polynomials.


```python
T = sp.Function('T')
Tminus1 = sp.Function('Tminus1')
Tminus2 = sp.Function('Tminus2')
```


```python
def Che(r):
    Tminus2 = 1
    Tminus1 = x
    if r == 0:
        return [Tminus2]
    elif r == 1:
        return [Tminus2,Tminus1]
    else:
        Tlist = [Tminus2,Tminus1]
        for i in range(r-1):
            T = 2*x*Tminus1 - Tminus2
            T = sp.simplify(T)
            Tlist.append(T)
            Tminus2 = Tminus1
            Tminus1 = T
    return Tlist
```

Then we verify that the polynomials are the same.


```python
Chelist = Che(3)
print(Chelist)
Chelistdivided = [1] + [sp.simplify(Chelist[i]/(2**(i - 1))) for i in range(1, len(Chelist))]
print(Chelistdivided)
tlist = ortho(1/sp.sqrt(1-x**2),-1,1,3)
print(tlist)
print(tlist == Chelistdivided)
```

    [1, x, 2*x**2 - 1, x*(4*x**2 - 3)]
    [1, x, x**2 - 1/2, x*(x**2 - 3/4)]
    [1, x, x**2 - 1/2, x*(x**2 - 3/4)]
    True


We may find the roots for the polynomials.


```python
sp.solve(tlist[2],x)
```




    [-sqrt(2)/2, sqrt(2)/2]



6d848739c053f4edcd94be584088cd2332a2de79
ipynb
Jupyter Notebook
M2AA3/M2AA3-Least-Square/Lesson 02 Best Approximation Problem/Orthogonal Polynomial.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
2021-05-16T04:08:34.000Z
2021-11-27T12:56:10.000Z
M2AA3/M2AA3-Least-Square/Lesson 02 Best Approximation Problem/Orthogonal Polynomial.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
M2AA3/M2AA3-Least-Square/Lesson 02 Best Approximation Problem/Orthogonal Polynomial.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
2020-03-31T00:23:13.000Z
2020-05-13T15:01:46.000Z
Qwen/Qwen-72B
1. YES
2. YES