Algorithms: Devising quadrature (integration) formulas

Very large scale projects such as modeling of atomic interactions with time-dependent fields, or FEM (finite element) analysis of electric and magnetic field geometries require computations of vast numbers of similar integrals with each iteration of the main relaxation loop. You can construct fast integration schemes given the types of functions that must be integrated that are specifically tailored to speed up your analysis.
Suppose that we we need to integrate linear functions only as a trivial example. We assume that for a first order polynomial, the formula

$\int_{x_)}^{x_1} \, f(x) \, dx= a f(x_0)+b f(x_1)$

is exact. Next compute the integral exactly

$\int_{x_)}^{x_1} \, f(x) \, dx= f(x_0)\big(x_1-x_0\big) +{1\over 2}\big(f(x_1)-f(x_0)\big)\big(x_1-x_0\big)$

$={1\over 2}\big(f(x_1)+f(x_0)\big)\big(x_1-x_0\big)$

and solve for the coefficients

$a=b={1\over 2}\big(x_1-x_0\big)={1\over 2}\delta$

and we arrive at

$\int_{x}^{x+\delta} \, f(x) \, dx={1\over 2}\delta f(x)+{1\over 2}\delta f(x+\delta)$

which we can iterate by gluing together many intervals to build up the domain $\{x \, | x_0\le x\le x_N \}$ from $N$ segments, with

$\delta={x_N-x_0\over N}, \qquad x_i=x_0+i\delta$

$\int_{x_0}^{x_N} \, f(x) \, dx=\sum_{i=1}^N{1\over 2}\delta \Big(f(x_{i-1})+f(x_i)\big)=\delta\big({1\over 2} f(x_0)+\sum_{i=1}^{N-1} f(x_i)+{1\over 2} f(x_N)\Big)$

This is a very good approximation if the function is linear on all of the little intervals

$[x_n,x_n+\delta]=[n\delta, (n+1)\delta]$

that we have chained together.

Next assume that $f(x)$ is quadratic on the interval $[x,x+2\delta]$ .

$\int_x^{x+2\delta} x^2 \, dx= a x^2 +b (x+\delta)^2 +c (x+2\delta)^2$

and construct the coefficients $a,b$ and $c$ from the resulting algebraic equations

${1\over 3} (x+2\delta)^3-{1\over 3} x^3=2\delta x^2+4\delta^2 x +{8\over 3}\delta^3=a x^2+b(x^2+2x\delta+\delta^2)+c(x^2+4\delta x +4\delta^2)$

This results in

$2\delta=a+b+c, \qquad 4\delta^2=2\delta b+4\delta c,\qquad {8\over 3}\delta^3=b\delta^2+4\delta^2 c$

$a={\delta\over 3}, \qquad b={4\delta\over 3}, \qquad c={\delta\over 3}$

We have obtained Simpson’s rule, which is better than trapezoid for larger intervals

$\int_x^{x+2\delta} f(x) \, dx={\delta\over 3}\Big(f(x)+4 f(x+\delta)+f(x+2\delta)\Big)$

which is exact for any second order $f(x)$ . Differentiate with respect to $x$ and we get the variant of Simpson’s rule that we used in the Adams-Bashford-Milne integration algorithm

$f(x+\delta)-f(x)={\delta\over 3}\Big(f'(x)+4 f'(x+\delta)+f'(x+2\delta)\Big)$

Rename:

$\int_x^{x+\delta} f(x) \, dx={\delta\over 6}\Big(f(x)+4 f(x+\delta/2)+f(x+\delta)\Big)$

and if $f(x)$ is quadratice on $[a,b]$ then

$\int_a^{b} f(x) \, dx={(b-a)\over 6}\Big(f(a)+4 f({a+b\over 2})+f(b)\Big)$

and you can chain together adjacent intervals upon which $f(x)$ satisfies this. You can see that the resulting formulas are very similar to Gaussian quadrature schemes in which we evaluate the integral with a weighted sum of the function values at selected points on the interval.

Integrals such as $\int_{-\infty}^\infty f(x) \, dx$ can be handled by mapping $[-\infty,\infty]\rightarrow [-1,1]$ and then using Gaussian quadrature http://abadonna.physics.wisc.edu/2021/08/01/algorithms-for-numerical-integration/ or Clenshaw-Curtis: let $x=L\cot(y)$ , for some well-chosen $L$ (in most cases $L=1$ works well)

$\int_{-\infty}^\infty f(x) \, dx=L\int_0^\pi {f(L\cot(y)) dy\over \sin^2y}\approx {L\pi\over N}\sum_{k=1}^{N-1} {f(L\cot(k\pi/N))\over \sin^2(k\pi/N)}$

In julia…

Clenshaw2sided=function(f,N,L)
abscissas=(k*π/N for k in 1:(N-1))
weights=sin.(abscissas)
weights=weights.^(-2)
args=f.(L*cot.(abscissas))
(L*π/N)*sum(weights .* args)
end

For Gaussian quadrature of $\int_{-1}^1 {f(x)\over\sqrt{1-x^2}}dx$ using Chebyshev polynomials:

quad=function(f,N)
abscissas=((2.0*k+1)*π/(4.0*N+2) for k in 0:N)
weights=abs.(sin.(abscissas))
args=f.(cos.(abscissas))
(π/(N+1.0))*sum(weights .* args)
end

Home 2.0

Algorithms: Devising quadrature (integration) formulas

Related Posts

Algorithms: Devising quadrature (integration) formulas

Published by Jeff on September 18, 2021September 18, 2021

Related Posts

Quantum spin-chains-III. Density Matrix Renormalization Group

Protected: Autodifferentiation with dual numbers

Protected: Number theory and combinatorics for classical encryption-I