Quadrature¶

Trapezoid integration¶

• $N_{trap}$ is the number of trapezoids ($N_{trap}=5$).
• We have $N_{trap}+1$ points (6 points).

\begin{align} I_0 &= \frac{\Delta x}{2}(f_0+f_1), \\ I_1 &= \frac{\Delta x}{2}(\phantom{f_0+}f_1+f_2), \\ I_2 &= \frac{\Delta x}{2}(\phantom{f_0+f_1+}f_2+f_3), \\ I_3 &= \frac{\Delta x}{2}(\phantom{f_0+f_1+f_2+}f_3+f_4), \\ I_4 &= \frac{\Delta x}{2}(\phantom{f_0+f_1+f_2+f_3+}f_4+f_5). \\ \end{align}$$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$:

\begin{align} 1:&\phantom{xxxxx}I_1 = \int_{-1}^11dz = \left.z\right|_{-1}^1 &= 2 &\Leftrightarrow \sum_{i=0}^1w_iz_i^0 = w_0z_0^0 + w_1z_1^0, \\ z:&\phantom{xxxxx}I_z = \int_{-1}^1zdz = \left.\frac{1}{2}z^2\right|_{-1}^1 &= 0 &\Leftrightarrow \sum_{i=0}^1w_iz_i^1 =w_0z_0 + w_1z_1, \\ z^2:&\phantom{xxxxx}I_{z^2} = \int_{-1}^1z^2dz = \left.\frac{1}{3}z^3\right|_{-1}^1 &= \frac{2}{3} &\Leftrightarrow \sum_{i=0}^1w_iz_i^2 =w_0z_0^2 + w_1z_1^2, \\ z^3:&\phantom{xxxxx}I_{z^3} = \int_{-1}^1z^3dz = \left.\frac{1}{4}z^4\right|_{-1}^1 &= 0 &\Leftrightarrow \sum_{i=0}^1w_iz_i^3 =w_0z_0^3 + w_1z_1^3. \\ \end{align}

• This gives four equations in four unknowns: $z_0$, $z_1$, $w_0$, $w_1$:

\begin{align} 2 &= w_0 + w_1 \\ 0 &= w_0z_0 + w_1z_1 \\ \frac{2}{3} &= w_0z_0^2 + w_1z_1^2 \\ 0 &= w_0z_0^3 + w_1z_1^3 \end{align}

• Solving these gives:

\begin{align} w_0 &= 1, \\ w_1 &= 1, \\ z_0 &= -\frac{1}{\sqrt{3}}, \\ z_1 &= \frac{1}{\sqrt{3}}. \end{align}

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

Numpy has a function that gives the weights and abscissas.¶

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.$$

Example¶

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

• 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

\begin{align} M_0 &= \int m^0n(m)dm = m_0^0w_0 + m_1^0w_1, \\ M_1 &= \int m^1n(m)dm = m_0^1w_0 + m_1^1w_1, \\ M_2 &= \int m^2n(m)dm = m_0^2w_0 + m_1^2w_1, \\ M_3 &= \int m^3n(m)dm = m_0^3w_0 + m_1^3w_1. \end{align}

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.