Approximate Numerical Convolution with a Singularity in the kernel

Question

Although similar to the other Convolution of function with singularity question, it is different as I am suggesting an approximation method but asking how to make it more rigorous. I asked on Signal Processing SE but it was suggested I ask here instead. The question there has more plots that I could not include here due to low reputation points.

I am trying to perform a numerical convolution where one of the operands is singular at t = 0. Concretely:

$\int_{0}^{t}G(t-\tau)\dot{\epsilon}(\tau)d\tau$,

where $\dot{\epsilon}(\tau)$ can be assumed a 'nice' function (smooth, no singularities etc.) but $G$ has a singularity at 0 and generally takes the form of something like

$G(t) = t^{-\beta}$, with $0 <\beta < 1$,

but could be a bit more complicated. I have solved this analytically for some 'nice' $\dot{\epsilon}(\tau)$, e.g. $\dot{\epsilon}(\tau)=\tau$ and although the integral is improper it converges. Ultimately $\dot{\epsilon}(\tau)$ will need to be experimental data so I need to solve the convolution numerically, and will not necessarily be able to change the way time is discretized. It cannot be solved using standard FFT methods due to the singularity and using quadrature is slow and also fails due to the singularity.

The current numerical approach I'm trying is just approximating the singular point by a well defined value, e.g. letting $G(0)=G(\Delta t/2)$. This works well sometimes, but the accuracy of approximation depends on $\beta$ and the time step $\Delta t$.

So my questions:

1. If the above approximation is OK, what is a more rigorous way of approximating the point at 0 that will give the best approxomation for any singular $G(t)$ with any time step and $\beta$?

2. If not, what is the fastest, correct way to compute the convolution with a singularity?

Any textbooks or papers as references would be very welcome, though the simpler and clearer the explanation the better.

Answer 1

Use of numerical quadrature for singular integrals is a fairly significant area of active research, as they can be used to discretize and thus solve integral equations that are used in modeling a variety of problems in physical science.

One general strategy is, if you know the asymptotics of the singularity at $x = 0$, to separate the integral into two pieces. Away from the singularity, you can use standard quadrature rules that are accurate for very smooth functions. Near the singularity, use the known asymptotics of the singularity (for example, if you know that the integrand grows like $|x|^{- \beta}$ as you describe) to form a new quadrature that takes advantage of the exact known integral for $|x|^{-\beta}$.

For example, consider the computation of $$ I = \int_0^1 x^{-1/2} f(x) \, dx, $$ where $f(x)$ is analytic. Then locally about $x = 0$, the integrand looks like $x^{-1/2} (f(0) + x f'(0) + O(|x|^2))$.

For small $\epsilon$, we use $$ \int_0^{\epsilon} x^{-1/2} f(x) \, dx = \int_0^{\epsilon} x^{-1/2} (f(x) - f(0)) \, dx + \int_0^{\epsilon} x^{-1/2} f(0) \, dx. $$ The first integrand behaves like $f'(0) x^{1/2} + O(\epsilon^{3/2})$ and can be computed using an appropriate quadrature for such an integrand, while the second integral is computed exactly.