Quadrature

Trapezoid integration

• 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:

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

• 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:

• (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:

• We transform $x$ to $z$ using:

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

• Hence,

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

Example:

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

• 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
• Gauss-Chebyshev: $$W(z)=(1-z^2)^{-1/2},$$ $$-1
• Gauss-Laguerre: $$W(z) = z^{\alpha}e^{-z},$$ $$-\infty
• Gauss-Hermite: $$W(z)=e^{-z^2},$$ $$-\infty
• Gauss-Jacobi: $$W(z) = (1-z)^{\alpha}(1+z)^{\beta},$$ $$-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

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

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

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.