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 Trapazoid rule, $w_i=\Delta x$ for the interior points and and $w=\Delta x/2$ for the edge points.
Gauss Quadrature selects the $x_i$ and the $w_i$ so that the integrals of polynomials up to degree $2n-1$ is exact.
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:
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: 1.6928356957111381 I Python quad func: 1.6928356957111381
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.]
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
def f(x): # define the function
#return np.exp(x)
#return np.sin(x)
#return x**8
#return x**(2-1/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 = quad(f,a,b)[0] # "exact" integral: quad function
n = np.arange(2,10) # 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 trapazoid 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(['Trapazoid', 'Gauss quadrature'], frameon=False);
$$-\infty<x<\infty$$
This is given by
$$S(n) = n(m)K_gm^{-2/3}.$$Here, $K_g$ is a constant. the source term for $M_1=\rho Y_s$ will have units of $kg/m^3s$. (The power here is really -1/3, but we'll use -2/3 for illustration) The source term becomes
$$\int m^kS(n)dm \rightarrow \int n(m)m^kK_gm^{-2/3}dm $$ $$\rightarrow \int n(m)K_gm^{k-2/3}dm.$$