Quadrature

Trapezoid integration

$$I = \left(\Delta x\sum_{i=0}^{5}f_i\right) -\frac{\Delta x}{2}(f_0+f_{5}).$$

$$I = \left(\Delta x\sum_{i=0}^{N_{trap}}f_i\right) -\frac{\Delta x}{2}(f_0+f_{N_{trap}}).$$

Note the overlap of the terms: two of every $f$ value except first and last.
- So, we double all of them and then subtract the first and last.

This has the following generic form, which is common to other numerical integration (quadrature) methods:

$$I = \int_a^bf(x)dx \approx \sum_{i=0}^{n-1}w_if(x_i).$$

That is, the integral is approximated as some sum over the function at points $x_i$, multiplied by some coefficients $w_i$. For the Trapezoid rule, $w_i=\Delta x$ for the interior points and and $w=\Delta x/2$ for the edge points.

Usually the points $x_i$ are fixed, and a uniform spacing is common.
- Call these $x_i$ points the abscissas
Here, we allow both specification of the weights $w_i$ and the abscissas $x_i$.
This will allow much better accuracy for a given computational cost.

Gauss Quadrature selects the $x_i$ and the $w_i$ so that the integrals of polynomials up to degree $2n-1$ is exact.

This is done on the domain $[-1,1]$ instead of $[a,b]$ for generality and convenience.
We will write $z$ as our integration variable instead of $x$ since the quadrature is written on the domain $[-1,1]$. We will then transform back to $[a,b]$ and use the more common symbol $x$.

Two point quadrature: $n=2$

Integrate polynomials to degree $2n-1=3$
Polynomials: $1$, $z$, $z^2$, $z^3$:
This gives four equations in four unknowns: $z_0$, $z_1$, $w_0$, $w_1$:
Solving these gives:
Hence, for some $p(z) = a_3z^3 + a_2z^2 + a_1z + a_0 dz$, for some $a$’s, we have

$$\int_{-1}^1 p(z)dz = 1\cdot p(-1/\sqrt{3}) + 1\cdot p(1/\sqrt{3}).$$

This can be extended to higher degree using higher $n$.
To the extent that our function is well represented by a degree $2n-1$ polynomial, the following approximation is accurate:

$$\int_{-1}^1 F(z) dz \approx \sum_0^{n-1} w_iF(z_i).$$

(Recall, a Taylor Series approximation of a function is a polynomial representation of the function.)

Example: demonstrate using a random 4th order polynomial

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
%matplotlib inline

a = np.random.rand(4)

Igq = np.polyval(a, -1.0/np.sqrt(3)) + np.polyval(a, 1.0/np.sqrt(3))
Ipy = quad(lambda z: np.polyval(a, z), -1, 1)[0]

print('I Gauss Quadrature: ', Igq)
print('I Python quad func: ', Igq)

I Gauss Quadrature:  2.292327052718544
I Python quad func:  2.292327052718544

Numpy has a function that gives the weights and abscissas.

from numpy.polynomial.legendre import leggauss

n    = 2
z, w = leggauss(n)

print('abscissas: z_i', z)
print('weights  : w_i', w)

abscissas: z_i [-0.57735027  0.57735027]
weights  : w_i [1. 1.]

Scale the domain

We often want to solve on some domain $[a,,b]$ instead of $[-1,,1]$.
That is, suppose we have a function $f(x)$ we want to integrate:

$$I = \int_a^bf(x)dx.$$

We transform $x$ to $z$ using:

$$x = mz+c.$$

At $x=a$, $z=-1$, and at $x=b$, $z=1$.
This gives

$$m=\frac{b-a}{2},$$ $$c = \frac{b+a}{2}.$$

Hence,

$$ x = \left(\frac{b-a}{2}\right)z + \frac{b+a}{2}.$$

Now change the integration variable from $x$ to $z$:

$$ I = \int_a^bf(x)dx = \int_{-1}^1f(mz+c)mdz = m\int_{-1}^1f(mx+c)dz.$$

$$I = \int_a^bf(x)dx = m\int_{-1}^1f(mz+c)dz,$$

Example:

$$I = \int_a^b f(x)dx = \int_3^7 4+[x-5]^3dx,,\text{or,}$$ $$a = 3,,,,b = 7,$$ $$x = mz+c,$$ $$m = \frac{b-a}{2} = \frac{7-3}{2} = 2,$$ $$c = \frac{b+a}{2} = \frac{7+3}{2} = 5,$$ $$I = m\int_{-1}^1 f(mz+c)dz = m\int_{-1}^1 4+[(mz+c)-5]^3dz.$$

def f(x):
    return 4+(x-5)**3

a = 3
b = 7
m = (b-a)/2
c = (b+a)/2


I_fxab =   quad(f,                   a, b)[0]
I_fz11 = m*quad(lambda z: f(m*z+c), -1, 1)[0]

print(f"Compare integral over a<x<b to integral over -1<z<1 using quad function: {I_fxab}, {I_fz11}")

#--------------

n = 2
z, w = leggauss(n)
I_GQ = m*np.sum(w*f(m*z+c))

print("Compare to Gauss Quadrature: ", I_GQ)

Compare integral over a<x<b to integral over -1<z<1 using quad function: 16.0, 16.0
Compare to Gauss Quadrature:  16.0

Example

Compare the error of integration for Gauss Quadrature and Trapezoid methods for varying number of quadrature points.

def f(x):                   # define the function
    #return np.exp(x)
    return np.sin(x)
    #return x**8
    #return x**(5/3)
    #return np.sqrt(x)
    #return x**(-1/3)

a = 3
b = 7
m = (b-a)/2
c = (b+a)/2
#Iexact = np.exp(b)   - np.exp(a)                # exact integral
Iexact = -np.cos(b)  + np.cos(a)                # exact integral
#Iexact = 1/9*b**9    - 1/9*a**9                   # exact integral
#Iexact = 3/8*b**(8/3)- 3/8*a**(8/3)        # exact integral
#Iexact = 2/3*b**(3/2)- 2/3*a**(3/2)        # exact integral
#Iexact = 3/2*b**(2/3) - 3/2*a**(2/3)       # exact integral
#Iexact = quad(f,a,b)[0]                       # "exact" integral: quad function

n = np.arange(2,18)             # list of number of quadrature points to try
ncases = len(n)                 # number of cases
Ierr = np.empty(ncases)         # array of relative errors

for i,ni in enumerate(n):
    z, w = leggauss(ni)     
    Ierr[i] = np.abs((m*np.sum(f(m*z+c)*w) - Iexact)/Iexact)
    
#-------------- Compare to trapezoid method with the same number of points:

def Itrap(f, a,b,n):
    x = np.linspace(a,b,n)
    dx = x[1] - x[0]
    return dx*np.sum(f(x)) - dx/2*(f(x[0])+f(x[-1]))

Ierr_trap = np.empty(ncases)
for i, ni in enumerate(n):
    Ierr_trap[i] = Itrap(f, a, b, ni)

Ierr_trap = np.abs( (Ierr_trap - Iexact)/Iexact )

#------------- plot
    
plt.rc('font', size=14)
plt.plot(n, np.log10(Ierr_trap), 'o-');
plt.plot(n, np.log10(Ierr), 'o-');
plt.xlabel("number of points");
plt.ylabel(r'$\log_{10}$ relative error');
plt.legend(['Trapezoid', 'Gauss quadrature'], frameon=False);

The linear convergence of Gauss Quadrature, when plotted on $\log(\epsilon_r)$ versus $n$ axes implies exponential convergence rate of Gauss Quadrature. Trapezoid is only a power law.

Generalized quadratures

Often, quadratures are written more generally as

$$\int_{-1}^1W(z)F(z)dz \approx\sum_{i=0}^{n-1}w_if(z_i).$$

Here, $W(z)$ is a weight function.
Above, we used $W(z)=1$. This is called a Gauss-Legendre Quadrature.
Guass-Legendre Quadrature is only accurate to the extent that the function can be well represented by a polynomial over the integration domain.
However, it may happen that our function can be split into a weight function multiplied by a polynomial-like part: $$F(z) \rightarrow W(z)\hat{F}(z).$$
As Numerical Recipes states: “We can arrange the choice of weights and abscissas to make the integral exact for a class of integrands polynomials times some known function $W(z)$ rather than for the usual class of integrands polynomials.”
When solving for the weights $w_i$ and abscissas $z_i$, rather than integrating the polynomials $1$, $z$, $z^2$, etc., as above, we integrate $W(z)1$, $W(z)z$, $W(z)z^2$, etc.
- This assumes we can do these integrations.

Quadrature types

Guass-Legendre: $$W(z)=1,$$ $$-1<x<1$$
Gauss-Chebyshev: $$W(z)=(1-z^2)^{-1/2},$$ $$-1<x<1$$
Gauss-Laguerre: $$W(z) = z^{\alpha}e^{-z},$$ $$-\infty<x<\infty$$
Gauss-Hermite: $$W(z)=e^{-z^2},$$ $$-\infty<x<\infty$$
Gauss-Jacobi: $$W(z) = (1-z)^{\alpha}(1+z)^{\beta},$$ $$-1<x<1$$

Example: particle size distributions

soot formation
aerosols
cloud droplets
sprays
griding

Consider a particle size distribution $n(m)$, where $m$ is mass of a particle, and $n$ is number of particles per unit volume per unit mass.

Then we call this a particle density function.
$N = \int n(m)dm$, where $N$ is the total number of particles per unit volume.

Approaches

Direct method: track all possible sizes $m$. But this is expensive, there may be millions of sizes.
Sectional method: track just course bins, say 100.
Moment methods, just track the first few moments of the distribution

Moment methods

The $k^{th}$ moment of the distribution is given by

$M_k = \int n(m)m^k dm$.

(Notice the resemblance to a weighted Gauss quadrature if $n(m)$ is our weight function.

We care most about the first two moments $M_0$ and $M_1$.

For $k=0$, $M_0=N$ is just the total particle number density (#/m$^3$).
For $k=1$, $M_1=\rho Y_p$ is just the particle mass density (kg of particles per m$^3$).

The more moments that are solved, the more accurate the results, especially $M_0$ and $M_1$.

Suppose we know how the size distribution $n(m)$ evolves. That is, we know

$$\frac{dn(m)}{dt} = S(m, n(m)).$$

The source term $S$ accounts for collisions, chemical growth, etc. (If two particles collide and stick, then $n(m)$ shifts to have more particles at larger sizes.)

We can convert this into an equation for $dM_k/dt$ by multiplying by $m^k$ and integrating over all sizes. This gives

$$\frac{dM_k}{dt} = \int m^kS(m, n(m))dm.$$

Problem

So, we evolve the moments in time according to this equation.

But, we need $n(m)$, and we need to be able to do the integral.
But, the whole point of using moments was to avoid knowing $n(m)$ and evolving it!
This is called a closure problem.

Solution

We treat the integral with quadrature, using $n(m)$ as our weight function.

The source integral is then

$$\int m^kS(m, n(m))dm \approx \sum_im_i^kS(m_i, w_i).$$

Now, we need the $m_i$ and the $w_i$.

But we can evaluate these from the moments $M_k$, which we are tracking.
Suppose we use two quadrature points. Then we have two weights $w_0$ and $w_1$, and two abscissas $m_0$ and $m_1$.
That’s four unknowns. We need four equations.
Track four moments $M_0$, $M_1$, $M_2$, $M_3$.

Using the definition of the moments we have

Solve this system for $m_0$, $m_1$, $w_0$, $w_1$, then we can compute the integrals in the $dM_k/dt$ equations and the system is closed.

Note how we evaluate the moment source terms in terms of $n(m)$ without ever having to know $n(m)$. This works because we are dealing with integrals using a weighted Gauss quadrature where we $do$ know the moments. Those moments give use the weights and abscissas for the quadrature, allowing us to do the integrations.

This is an extremely accurate and powerful method.