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:
using Plots
using Random
using QuadGK
using Polynomials
using CurveFit
using FastGaussQuadrature
using PyFormattedStrings
using LaTeXStrings
a = rand(4)
Igq = Polynomial(a)(-1.0/sqrt(3)) + Polynomial(a)(1.0/sqrt(3))
Ipy = quadgk(z -> Polynomial(a)(z), -1, 1)[1]
println("I Gauss Quadrature: $Igq")
println("I Python quad func: $Ipy")
I Gauss Quadrature: 1.5055565082103142 I Python quad func: 1.5055565082103142
n = 2
z, w = gausslegendre( n );
println("abscissas: z_i $z")
println("weights : w_i $w")
abscissas: z_i [-0.5773502691896258, 0.5773502691896258] weights : w_i [1.0, 1.0]
f(x) = 4+(x-5)^3
a = 3
b = 7
m = (b-a)/2
c = (b+a)/2
I_fxab = quadgk(f, a, b)[1]
I_fz11 = m*quadgk(z -> f(m*z+c), -1, 1)[1]
println(f"Compare integral over a<x<b to integral over -1<z<1 using quad function: {I_fxab}, {I_fz11}")
#--------------
n = 2
z, w = gausslegendre(n)
I_GQ = m*sum(w.*f.(m*z.+c))
println("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
function f(x)
#return exp(x)
return sin(x)
#return x^8
#return x^(5/3)
#return sqrt(x)
#return x^(-1/3)
end
a = 3
b = 7
m = (b-a)/2
c = (b+a)/2
#Iexact = exp(b) - exp(a) # exact integral
Iexact = -cos(b) + 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 = quadgk(f,a,b)[1] # "exact" integral: quad function
n = 2:16 # list of number of quadrature points to try
ncases = length(n) # number of cases
Ierr = zeros(ncases) # array of relative errors
for (i,ni) in enumerate(n)
z, w = gausslegendre(ni)
Ierr[i] = abs((m*sum(f.(m*z.+c).*w) .- Iexact)/Iexact)
end
#-------------- Compare to trapazoid method with the same number of points:
function Itrap(f, a,b,n)
x = LinRange(a,b,n)
dx = x[2] - x[1]
return dx*sum(f.(x)) .- dx/2*(f(x[1])+f(x[end]))
end
Ierr_trap = zeros(ncases)
for (i, ni) in enumerate(n)
Ierr_trap[i] = Itrap(f, a, b, ni)
end
Ierr_trap = abs.( (Ierr_trap .- Iexact)/Iexact );
plot( n, log10.(Ierr_trap), marker=:circle, lw=2, label="trapazoid")
plot!(n, log10.(Ierr), marker=:circle, lw=2, label="Gauss quadrature")
plot!(xlabel="# points")
plot!(ylabel=L"$\log_{10}$ relative error")
plot!(legend_foreground_color = nothing)
KERNEL EXCEPTION KeyError: key "debug_request" not found Stacktrace: [1] getindex(h::Dict{String, Function}, key::String) @ Base ./dict.jl:498 [2] eventloop(socket::ZMQ.Socket) @ IJulia ~/.julia/packages/IJulia/6TIq1/src/eventloop.jl:8 [3] (::IJulia.var"#14#17")() @ IJulia ./task.jl:484
$$-\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.$$