I have a set of PDEs governing some function $f(r)$ which I desire to solve via a psuedospectral method (we can consider $f$ to be smooth). It is defined on the interval $r\in[0,\infty)$ with symmetry at $r=0$ and $f\to f_\infty$ (a constant) as $r\to\infty$. Accordingly, I suspect something like a rational Chebyshev basis would be appropriate here.
Regardless of the basis functions chosen, I will be able to determine $f(r)$ at a set of a grid points $\{r_i\}_{i=1}^N$. Part of the governing PDEs is an integral of the form $$I(r)=\int_0^\infty K(r',r)f'(r')\,{dr}'$$ which I need to know the value of at each grid node. The kernel $K$ is symmetric and singular at $r'=r$, so I split the integral into two parts: $$I(r)=\int_{0}^r K_1(r',r)f'(r')\, {dr}'+\int_{r}^\infty K_1(r,r')f'(r')\, {dr}'.$$ $K_1$ is also singular at $r'=r$ with asymptotic form $K_1(r',r)\sim\log(1-(r'/r)^2)$ (logarithmic singularity).
My question is how I should evaluate these integrals via a quadrature that retains the spectral accuracy of the numerical method. If the kernels were not singular, this would be very simple with the quadrature points and weights for the set of basis functions chosen. Essentially, I wish to write $$I(r_j)=\sum_{i=1}^{j}K_1(r_i,r_j)f'(r_i)w_i+\sum_{i=j}^{N}K_1(r_j,r_i)f'(r_i)w_i.$$ Now, it is not this simple since the integrand is singular when $i=j$. I have seen a few methods used in literature:
- Using a tanh mapping $r'/r=\tanh(z/L)$ transforms the interval $r'\in(0,r)$ to $z\in(0,\infty)$ and the integrand becomes nonsingular. It is then recommended to use a spectrally accurate quadrature on this new integral of the form $$I(r_j)=\frac{r_j}{L}\int_0^\infty K_1(r_j\tanh(z/L),r_j)f'(r_j\tanh(z/L))\text{sech}^2(z/L)\,dz.$$ The hyperbolic secant here is sufficient to cause decay of the integrand as $z\to\infty$. The issue here is that the nodes at which I have the value of $f'$ are now $\{L\,\text{atanh}(r_i/r_j)\}_{i=1}^j$, which not only includes $\infty$, but also does not correspond to any spectral quadrature that I know of. The second integral can be transformed in a similar fashion with similar issues. Is there a way a mapping like this can be used successfully? Do I need to use computationally-expensive off-grid interpolation for this?
- Singular basis functions can be added to the basis set to capture the singularity better. However, I'm having a hard time finding any digestible references on this since it seems to be very problem specific. Is this a viable method to approaching this integral?
- Perform the integral analytically? If the basis functions for $f'$ are known, it might be possible to analytically perform the integral against the known logarithmic form of the singularity. I highly doubt such a nice analytical result exists, however. I have been unable to find one.
I would appreciate any insight you can give me!
The best way I have found to do this is as follows. If the function $f$ can be expressed as an expansion of the form $$f(r)=\sum_i a_i\phi_i(r)$$ for some basis set $\{\phi_i\}$ (this could be the cardinal functions using the known set of values), then the integral of interest becomes $$I(r)=\int_0^\infty K(r,r')f(r')dr'=\sum_ia_i\int_0^\infty K(r,r')\phi_i(r')dr'=\sum_ia_im_i(r).$$ where we denote the "mass" $$m_i(r)=\int_0^\infty K(r,r')\phi_i(r')dr'.$$ Since the functional forms of both $K$ and $\phi_i$ are known, these masses can be computed at each grid point $r_j$ to machine-precision by adaptive quadrature and then the sum for $I(r_j)$ evaluated. The convergence of the approximation of this integral is then the same as as that of the approximation of $f$ (which should be spectrally accurate). (I had an issue where the mass integrals did not converge as written above, so I had to shuffle around a factor of $r^2$ to make this work.)