In many quadrature problems, we are interested in computing $\int_a^b f(x) \, \mathrm{d}x$ via a quadrature sum. However, most software packages precompute the quadrature nodes and weights for use with integrals on $[-1,1]$ (Gaussian, Chebyshev and tanh-sinh quadratures are examples). Hence it is common to use a variable transformation $$ \int_{a}^{b} f(x) \, \mathrm{d}x = \frac{(b-a)}{2} \int_{-1}^{1} f(((b-a)x+(b+a))/2) \, \mathrm{d}x $$ The quadrature nodes are mapped from $x_{i} \in (-1,1) \mapsto t_{i} := ((b-a)x+(b+a))/2$. However, in floating point arithmetic, we can have $x_0 > -1$, but $t_0$ bitwise equal to $a$. If $f$ is singular at $a$, then this is catastrophic.

Is it possible to give a useful characterization of when interior nodes in the $(-1,1)$ domain will be squashed into boundary in the $[a,b]$ domain using this transformation?

Is there a well-known method of mitigating this problem?

Note: If $f$ is singular at $a$, I've been using $$ \int_{\mathrm{floatafter}(a)}^{b} f \, \mathrm{d}x, $$ but certainly this is not always necessary.

Edit: The goal here is for a quadrature routine to do the best with what it is given. I also assume the quadrature is adaptive.


1 Answer 1


In short:

Usually, the vanilla quadrature-based integration is used for well-behaved functions that do not have singularities on the integration domain. While some functions could sometimes still be integrated with that method relatively well, the singularity requires special treatment. When you are integrating a singular function, having floating point error for $x_0$ and $a$ in your notation is the least of your problems.

So, if you have a singularity at $a$, but you have to directly evaluate the function at the "floating point distance" close to it – you already have a problem. There is a lot of literature devoted to singularity treatment of different types for various applications.

Now, regarding the singularity treatment in more detail. The usual suspects:

  1. Extract the singularity and split the integral into two parts: singular and non-singular. The non-singular part can be integrated to the desired precision using standard quadratures, and non-singular part evaluated analytically.
  2. Alternatively, one can use specially developed quadrature rules for particular singular integrals. The nicest reference I have and commonly use is J. H. Ma, V. Rokhlin, and S. Wandzura, "Generalized Gaussian quadrature rules for systems of arbitrary functions," SIAM J. Numer. Anal., vol. 33, no. 3, pp. 971–996, Jun. 1996.
  3. Adaptive integration. The brute-force approach is usually the easiest to implement and often gives decent accuracy. Unfortunately, for highly singular functions it might never converge and is commonly way slower compared to the correctly formulated and implemented previous approaches.

Since I am from computational electromagnetics community, my examples and references are also from this area that is quite rich in singularities.

Example of the first approach in 2-D (common for computational electromagnetics):

$$ \int\limits_{C}H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)dl^\prime \tag{1} $$ In (1), we want to integrate the Hankel function of zero order and second kind over line segment $C$, $k$ is the wavenumber (for the purpose of this discussion can be considered as some constant), $\vec{\rho}$ is the constant position-vector to some observation point and $\vec{\rho}^\prime$ is the position vector on the line $C$ varying with $dl\prime$.

Now, if $\vec{\rho}$ is far away from the contour $C$, Hankel function behaves nicely; however, if $|\vec{\rho}-\vec{\rho}^\prime|\to0$ we have a singularity. The common treatment for this problem is to apply the small argument expansion (2) to $H_0^{(2)}$:

$$ H_0^{(2)}(k\underbrace{|\vec{\rho}-\vec{\rho}^\prime|}_{\to0})\to1-\frac{2j}{\pi}\left(\ln\frac{k|\vec{\rho}-\vec{\rho}^\prime|}{2}+\gamma\right) \tag{2} $$ where $j$ is the imaginary unity and $\gamma$ is Euler's constant. Now we can subtract expansion (2) from (1) and arrive to an expression that is ready for the numerical evaluation even if the observation point $\vec{\rho}$ lies on the integration contour $C$. $$ \int\limits_{C}H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)dl^\prime=\int\limits_{C}\underbrace{\left[H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)+\frac{2j}{\pi}\ln|\vec{\rho}-\vec{\rho}^\prime|\right]}_\text{*, non-singular}dl^\prime-\frac{2j}{\pi}\underbrace{\int\limits_C\ln|\vec{\rho}-\vec{\rho}^\prime|dl^\prime}_\text{**, singular} \tag{3} $$

In (3), the non-singular part (*) can be evaluated using regular quadrature, while for the singular logarithmic part the integral is evaluated analytically (the expression for analytic log-integration as well as the presented example is available from D. Wilton et al., "Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains," IEEE Trans. Antennas Propag., vol. 32, no. 3, pp. 276–281, Mar. 1984.).

Now, in many cases integration of (1) is possible using standard quadratures even for observation points $\vec{\rho}$ lying exactly on the line segment $C$. But one would certainly be required to use an expensive quadrature and suffer for the convergence.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.