# Gaussian Quadrature

## The idea behind Gaussian Quadrature
So far we focused on uniformly spaced grid points (nodes). Gaussian Quadrature is an algorithm
that chooses the nodes such that in the approximation
$$\int_a^b f(x)dx =\sum_{i=1}^{n}c_i f(x_i) \label{GQ}\tag{1}$$
the error is minimized in the sense that the result is _exact_ for the class of polynomials of degree $2n-1$.
This is because Gaussian Quadrature treats both $c_i$ and $x_i$ as parameters, i.e., there are $2n$ parameters which
is the number of parameters in a polynomial of degree $2n-1$. 

### Example: n=2 Gaussian Quadrature
Let's consider the case $n=2$ for $x\in [-1,1]$. We wish to determine $c_1, c_2, x_1, x_2$ such that
\begin{align}
\int_{-1}^{1} f(x)dx = &\sum_{i=1}^{n}c_i f(x_i) = c_1 f(x_1) + c_2 f(x_2) \label{condition}\tag{2}
\end{align}
is exact  whenever $f(x)$ is a polynomial of degree $2n-1=2\cdot 2-1=3$ or less, i.e.,
$$
f(x) = a_0 + a_1 x+a_2 x^2 + a_3 x^3.
$$

To achieve this we will first integrate this polynomial

$$
\int_{-1}^{1}f(x)dx = a_0 \int_{-1}^{1}dx + a_1 \int_{-1}^{1}x dx+ a_2 \int_{-1}^{1}x^2dx + a_3 \int_{-1}^{1}x^3dx.$$  

We can now consider the function $f(x)$ as a linear combination of the
4 functions $f_0(x)=1$, $f_1(x)=x$, $f_2(x)=x^2$, $f_3(x)=x^3$, i.e.,
$f(x)=a_0 f_0(x)+a_1 f_1(x)+a_2 f_2(x)+a_3 f_3(x)=\sum_{i=0}^3
a_if(x_i)$.  Thus, if for each $f_i(x),\ i=0,1,2,3$ we satisfy (\ref{condition}) because
\begin{split}
  \int_{-1}^{1}f(x)dx = & \int_{-1}^{1}\sum_i a_i f_i(x)dx =
  \sum_i a_i \int_{-1}^{1} f_i(x)dx = \sum_i a_i (c_1 f_i(x_1) + c_2
  f_i(x_2)) = \sum_i a_i c_1 f_i(x_1) + \sum_i a_i c_2 f_i(x_2)\\
   = &\  c_1 \sum_i a_i f_i(x_1) + c_2\sum_i a_i f_i(x_2) 
  = c_1 f(x_1) + c_2
  f(x_2).
\end{split}

Thus, for the 4 $f_i(x)$ polynomials we need to satisfy $\int_{-1}^1
f_i(x) dx = c_1 f_i(x_1) + c_2 f_i(x_2)$, which will provide 4
equations for the 4 parameters $c_1,c_2$, $x_1,x_2$:
  \begin{split}
    c_1 f_0(x_1) + c_2 f_0(x_2) = & c_1 \cdot 1 + c_2 \cdot 1 = \int_{-1}^1 1 dx = 2\\
    c_1 f_1(x_1) + c_2 f_1(x_2) = & c_1 \cdot x_1 + c_2 \cdot x_2 = \int_{-1}^1 x dx = 0\\
    c_1 f_2(x_1) + c_2 f_2(x_2) = & c_1 \cdot x_1^2 + c_2 \cdot x_2^2 = \int_{-1}^1 x^2 dx = \frac{2}{3}\\
    c_1 f_3(x_1) + c_2 f_3(x_2) = & c_1 \cdot x_1^3 + c_2 \cdot x_2^3 = \int_{-1}^1 x^3 dx = 0.
  \end{split}


This system can be solved easily by eliminating $c_1$ and $c_2$ from the first two equations, e.g.,
solve the first w.r.t $c_1$ and then replace the result in the second equation to express $c_2$ in
terms of $x_1$ and $x_2$. Then substitute $c_1$ and $c_2$ in terms of $x_1$ and $x_2$ in the last
two equations that will provide a system of 2 equations for the 2 unknowns $x_1$ and $x_2$. 

The solution
of the system is $c_1=1$, $c_2=1$, $x_1=-\sqrt{3}/3=-x_2$.  

**Therefore, the $n=2$ Gaussian quadrature
is
$$
\int_{-1}^{1} f(x)dx =  f(-\frac{\sqrt{3}}{3}) + f(\frac{\sqrt{3}}{3}),
$$
which yields an _exact_ result for any $f(x)$ that is a polynomial up to degree $n=3$!**



### <span style="color:blue"> Hands-on example: n=2 Gaussian quadrature on the [-1,1] interval </span>


```python
import numpy as np
#simple version of Gaussian Quadrature with n = 2
def GaussianQ_n2_unit_interval(f,*P):
    n = 2
    #define nodes (derived above)
    x = np.array([-np.sqrt(3)/3.,np.sqrt(3)/3.])
    #define coefficients (derived above)
    c = np.array([1.,1.])
    I = 0
    #now sum up the integral based on the n=2 Gaussian Quadrature formula
    for j in range(n):
        I += c[j]*f(x[j],*P)
    return I
```


```python
#define a fifth-order polynomial
def p_n5(x,a0,a1,a2,a3,a4,a5):
    return a0+a1*x+a2*x**2+a3*x**3+a4*x**4+a5*x**5

#integral of fifth-order polynomial over arbitrary integral
def int_p_n5(xmin,xmax,a0,a1,a2,a3,a4,a5):
    I = a0*(xmax-xmin) + a1*(xmax**2-xmin**2)/2. + a2*(xmax**3-xmin**3)/3. + a3*(xmax**4-xmin**4)/4. + a4*(xmax**5-xmin**5)/5. +a5*(xmax**6-xmin**6)/6.
    return I
#wrapper function for calling integral over [-1,1] interval
def int_p_n5_unit_interval(*P):
    return int_p_n5(-1,1,*P)
```


```python
#now test this by defining the coefficients of p_n5, our integrand
P = (1,1,1,0,0,0)
GQ_n2 = GaussianQ_n2_unit_interval(p_n5,*P)
exact = int_p_n5_unit_interval(*P)

print("exact integral: %e"%(exact))
print("Gaussian Quadrature: %e"%(GQ_n2))
```

    exact integral: 2.666667e+00
    Gaussian Quadrature: 2.666667e+00


<span style="color:blue"> Convince yourself that this is exact for `p_n5` with degree n<=3 </span>

### Arbitrary intervals
The general procedure for Gaussian quadrature outlined earlier applies
to the interval $[-1,1]$. However, we will often have general
integration domains $[a,b]$. We can still apply $\int_a^b f(x)dx =\sum_{i=1}^{n}c_i f(x_i)$ by a simple
transformation of the integration variable. An integral of the form
$$
\int_a^b f(x) dx 
$$
can easily be transformed to the form
$$\int_{-1}^1 f(t) dt$$

if we let

$$
t=\frac{2 x -a -b}{b-a}, 
$$

which implies
$$x=\frac{1}{2}[(b-a)t+a+b].
$$
Thus, when $x=a$, $t=-1$, and when $x=b$, $t=1$. We then have $dx=\frac{b-a}{2}dt$, and hence
$$
\int_a^b f(x) dx = \frac{b-a}{2}\int_{-1}^1 f\left(\frac{(b-a)t+(b+a)}{2} \right)dt.
$$
The last expression allows the application of our Gaussian Quadrature expression to
any function in any integration interval.


<span style="color:blue"> Now rewrite `GaussianQ_n2_unit_interval` to perform Gaussian Quadrature over any integration interval:
- write a function `tranform_variable(t, a,b)`
- use `tranform_variable(t, a,b)` within the n = 2 Gaussian Quadrature formula
- remember to apply the rescaling (b-a)/2 to the result obtained on the interval [-1,1]
</span>


```python
def transform_variable(t, a,b):
    return ((b-a)*t+a+b)/2
```


```python
def GaussianQ_n2_any_interval(f,a,b,*P):
    n = 2
    #define nodes on [-1,1] interval(derived above)
    t = np.array([-np.sqrt(3)/3.,np.sqrt(3)/3.])
    #define coefficients (derived above)
    c = np.array([1.,1.])
    I = 0
    #now sum up the integral based on the n=2 Gaussian Quadrature formula
    #evaluate f(x(t[j])) with x on the interval [a,b]
    for j in range(n):
        I += (b-a)/2.*c[j]*f(transform_variable(t[j],a,b),*P)
    return I
```


```python
P = (1,0,1,1,0,0)
a = 0
b = 2
GQ_n2 = GaussianQ_n2_any_interval(p_n5,a,b,*P)
exact = int_p_n5(a,b,*P)

print("exact integral: %e"%(exact))
print("Gaussian Quadrature: %e"%(GQ_n2))
```

    exact integral: 8.666667e+00
    Gaussian Quadrature: 8.666667e+00


### Gaussian Quadrature with n>2 
The same approach can yield results for degree $2n-1$ polynomials, but
this turns out to be more involved.  By using Legendre polynomials,
the process can be expedited. So, let's introduce the Legendre polynomials first.

**Legendre polynomials** ($P_n(x)$) are orthogonal ($\int_{-1}^1 P_n(x) P_m(x) = \delta_{nm}$) and form a complete basis.

Some of theirproperties are:
-  For each $n$, $P_n(x)$ is a monic polynomial of degree $n$, i.e., the coeficient of the smallest power of $x$
  that appears in the polynomial is $1$.
  
- $\int_{-1}^1 p(x) P_n(x) = 0$, when $p(x)$ is a polynomial of degree $<n$.
  
- The first 5 Legendre polynomials are: $P_0(x)$=1, $P_1(x)$=x, $P_2(x)=x^2-1/3$, $P_3(x)=x^3-3x/5$, \\ $P_4(x)=x^4-6x^2/7 +3/35$.

- <u>**Theorem**</u>
  Assume $x_1, x_2, \ldots x_n$ are the roots of the $n$-th Legendre
  polynomial $P_n(x)$ and for each $i=1,2,\ldots n$, consider the
  numbers
  $$
    c_i = \int_{-1}^{1} \prod_{j=1,j\neq i}^n \frac{x-x_j}{x_i-x_j} dx \label{eqci}\tag{3}
  $$
  If $p(x)$ is _any_ polynomial of degree less than $2n-1$, then
  $$
    \int_{-1}^1 p(x) dx = \sum_{i=1}^n c_i p(x_i).
  $$
 
 
 ### Example: n=3 Gaussian Quadrature
Let's consider the integral
$$
  I=\int_{-1}^1 e^x cos(x) dx = 1.933421496200713,
$$
and lets approximate it using Gaussian quadrature
with $n=3$, i.e.,

$$\tilde I = c_1 f(x_1) + c_2 f(x_2) + c_3 f(x_3).$$

where $x_1, x_2, x_3$ are the roots of the Legendre polynomial
$P_3(x)$, i.e., $x^3 - 3x/5 =0 \Rightarrow$ $x=0$ or $x^2-3/5=0
\Rightarrow$ $x=\pm \sqrt{3/5}$. Thus, $x_1=-\sqrt{3/5}$, $x_2=0$, and
$x_3=\sqrt{3/5}$. Next, we calculate the coefficients $c_i$ using Eq.(\ref{eqci}):
  \begin{split}
    c_1 = & \int_{-1}^{1} \prod_{j=1,j\neq 1}^n \frac{x-x_j}{x_1-x_j} dx = \int_{-1}^{1} \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3} dx = \frac{5}{9} \\
    c_2 = & \int_{-1}^{1} \prod_{j=1,j\neq 2}^n \frac{x-x_j}{x_2-x_j} dx = \int_{-1}^{1} \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3} dx = \frac{8}{9} \\
    c_3 = & \int_{-1}^{1} \prod_{j=1,j\neq 3}^n \frac{x-x_j}{x_3-x_j} dx = \int_{-1}^{1} \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2} dx = \frac{5}{9} \\
  \end{split}

Thus,
$$
\tilde I = \frac{5}{9} f\left(-\sqrt{\frac{3}{5}}\right) + \frac{8}{9} f(0) + \frac{5}{9} f\left(\sqrt{\frac{3}{5}}\right) = 1.933390469264298. \label{Gaussn4}\tag{4}
$$

Thus, the relative error is $|(\tilde I-I)/I|=1.6\times 10^{-5}$!!!  
**Notice that Eq.(\ref{Gaussn4}) applies to any function, and is _exact_ for any polynomial of degree 5 or less.**


### <span style="color:blue"> Hands-on example: n=3 Gaussian quadrature of $
I = \int_0^\pi x\cos(x) d x$</span>
- following the n = 2 example above, implement the n = 3 Gaussian quadrature expression on an arbitraty interval
- calculate the relative error by comparing to the exact value $I = 2$.


```python
def GaussianQ_n3_any_interval(f,a,b,*P):
    n = 3
    #define nodes on [-1,1] interval(derived above)
    t = np.array([-np.sqrt(3/5.),0,np.sqrt(3/5.)])
    #define coefficients (derived above)
    c = np.array([5./9,8./9,5./9])
    I = 0
    #now sum up the integral based on the n=3 Gaussian Quadrature formula
    for j in range(n):
        I += (b-a)/2.*c[j]*f(transform_variable(t[j],a,b),*P)
    return I
```


```python
def func(x,*P):
    return x*np.cos(x)
```


```python
exact_value = -2.
GQ_n3_approximation = GaussianQ_n3_any_interval(func,0,np.pi,*P)
fractional_error = np.abs((GQ_n3_approximation-exact_value)/exact_value)
print(fractional_error )
```

    0.0040611225242603854



```python

```

6b4be363b0ab574781371c4140039faf20d48c3a
ipynb
Jupyter Notebook
Lectures/Lecture 9/GaussianQuadrature_solutions.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
2020-09-10T06:45:46.000Z
2020-10-20T13:50:11.000Z
Lectures/Lecture 9/GaussianQuadrature_solutions.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Lectures/Lecture 9/GaussianQuadrature_solutions.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Qwen/Qwen-72B
1. YES
2. YES