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

So far we have seen Trapezoidal Rule and Simpson Rule. In general, we may approximate the integral $$I = \int_a^b w(x) f(x) dx$$
where $w \in C^0([a,b]), w>0$ and $w \in L^1((a,b))$ is a weight function, and $f \in L^1((a,b),w)$, by approximating $f$ as a polynomial. 

We may approximate $f$ with $n$-th order Largrange polynomial (sampling points $(x_i)_{i=0}^{k+1}$):
\begin{equation}
f(x) \approx \sum_{i=0}^n f(x_i) \ell_i(x) \quad \ell_i(x) = \prod_{0\leq j \leq i, \, j \neq i} \frac{x-x_j}{x_i-x_j}
\end{equation}

Then
\begin{equation}
\int_a^b w(x)f(x)dx \approx I_n = \sum_{i=0}^n f(x_i) \left(\int_a^b w(x) \ell_i(x) dx \right) = \sum_{i=0}^n w_i f(x_i) 
\end{equation}

We want to choose $(x_i)_{i=0}^n$ such that $I_n(f) = I(f)$ for $f$ polynomial of degree at most $m$, $m>n$ is maximised.

Turns out that if we choose roots of $\phi_{n+1}(x)$ (with $n+1$ degree), so that it belongs to the orthogonal basis of $L^2((a,b))$ (w.r.t. $\langle f,g \rangle = \int_a^b w(x) f(x)g(x) dx$, then $I_n(f) = I(f)$, $f$ any polynomial of degree at most $2n+1$. <br> 
Moreover, we can check by substituting $f = (\ell_i)^2$ that $w_i > 0$. We can find $w_i$ by such substitution, but this is usually tedious. Instead, we find $w_i$ by substituting $f = 1, x, x^2, ..., x_n$ and solve the simultaneous equation.

# Example 1
Let $f \in L^2([-1,1])$ and $w(x) = 1+4|x|$. Construct $I_2(f)$. 

1. Construct Orthogonal Basis: Recall the following formula:
$$\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$.


```python
# Define Symbols for inner product and norm
x = sp.symbols('x')
w = sp.Function('w')
f = sp.Function('f')
g = sp.Function('g')
```


```python
# inner product and norms
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))
```


```python
# Define symbols for basis
phi = sp.Function('phi')
phiminus1 = sp.Function('phiminus1')
phiminus2 = sp.Function('phiminus2')
```


```python
# Iteration
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
```


```python
# Compute the basis, up to order 2+1 = 3
olist = ortho(1+4*sp.Abs(x),-1,1,3) 
olist
```




    [1, x, x**2 - 4/9, x*(x**2 - 13/20)]



2. Compute the roots of $\phi_3(x)$, as $x_0, x_1, x_2$.


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




    [0, -sqrt(65)/10, sqrt(65)/10]



3. Compute $b_0 = I(1), \, b_1 = I(x), \, b_2 = I(x^2)$


```python
w = 1+4*sp.Abs(x)
b0 = sp.integrate(w,(x,-1,1))
b1 = sp.integrate(x*w,(x,-1,1))
b2 = sp.integrate((x**2)*w,(x,-1,1))
b = sp.Matrix([b0,b1,b2])
b
```




$\displaystyle \left[\begin{matrix}6\\0\\\frac{8}{3}\end{matrix}\right]$



4. Compute $w_0, w_1, w_2$ by solving
\begin{equation}
\begin{cases}
w_0 + w_1 + w_2 = b_0 \\
w_0 x_0 + w_1 x_1 + w_2 x_2 = b_1 \\
w_0 x_0^2 + w_1 x_1^2 + w_2 x_2^2 = b_2
\end{cases}
\end{equation}


```python
A = sp.Matrix([[1,1,1],[rootlist[0],rootlist[1],rootlist[2]], [rootlist[0]**2,rootlist[1]**2,rootlist[2]**2]])
A
```




$\displaystyle \left[\begin{matrix}1 & 1 & 1\\0 & - \frac{\sqrt{65}}{10} & \frac{\sqrt{65}}{10}\\0 & \frac{13}{20} & \frac{13}{20}\end{matrix}\right]$




```python
w0, w1, w2 = sp.symbols('w0,w1,w2')
sp.linsolve((A,b),[w0,w1,w2])
```




$\displaystyle \left\{\left( \frac{74}{39}, \  \frac{80}{39}, \  \frac{80}{39}\right)\right\}$



# Please edit below


```python
w = sp.exp(-x)
a = 0
b = sp.oo
n = 2
```


```python
# Compute the basis, up to order n+1
olist = ortho(w,a,b,3) 
olist
```




    [1, x - 1, (x - 3)*(x - 1) - 1, x**3 - 9*x**2 + 18*x - 6]




```python
# Compute root of phi(n+1)
rootlist = sp.solve(olist[2],x)
rootlist
```




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




```python
# Compute actual value of I(f) for f = 1, x,...
b0 = 1
b1 = 1
b = sp.Matrix([b0,b1])
b
```




$\displaystyle \left[\begin{matrix}1\\1\end{matrix}\right]$




```python
# Compute weight by solving equation
A = sp.Matrix([[1,1],[rootlist[0],rootlist[1]]])
w0, w1, w2 = sp.symbols('w0,w1,w2')
sp.linsolve((A,b),[w0,w1])
```




$\displaystyle \left\{\left( \frac{\sqrt{2}}{4} + \frac{1}{2}, \  \frac{1}{2} - \frac{\sqrt{2}}{4}\right)\right\}$




```python
# Integrate
integral = np.dot(np.array[w0], np.array[f(x) for x in rootlist])
```

08d67723466c3b54bd75812854c5eabc971f8554
ipynb
Jupyter Notebook
M2AA3/M2AA3-Polynomials/Lesson 04 - Quadrature/.ipynb_checkpoints/Gaussian Quadrature-checkpoint.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
2021-05-16T04:08:34.000Z
2021-11-27T12:56:10.000Z
M2AA3/M2AA3-Polynomials/Lesson 04 - Quadrature/.ipynb_checkpoints/Gaussian Quadrature-checkpoint.ipynb
ImperialCollegeLondon/Random-Stuff
219bc0e26ea6f5ee7548009c849959b268f54821
MIT
M2AA3/M2AA3-Polynomials/Lesson 04 - Quadrature/.ipynb_checkpoints/Gaussian Quadrature-checkpoint.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