Clenshaw-Curtis rules for oscillatory integrals with logarithmic singularities

This home web page for the project Clenshaw-Curtis rules for oscillatory with (possible) logarithmic singularities.

First, the papers where these rules were dealt:

• Domínguez, V. Filon-Clenshaw-Curtis rules for a class of highly-oscillatory integrals with logarithmic singularities. Journal of Computational and Applied Mathematics , 261 (2014), pp 299–319. Preprint available in http://arxiv.org/abs/1305.1365
• Domínguez, V, Graham, I.G., and Kim, T. Filon-Clenshaw-Curtis rules for highly-oscillatory integrals with algebraic singularities and stationary points. SIAM Journal of Numerical Analysis 51 (2013) no. 3, 1542–1566. Preprint available in http://arxiv.org/abs/1207.2283
• Domínguez, V., Graham, I., and Smyshlyaev, V. Stability and error estimates for Filon-Clenshaw-Curtis rules for highly-oscillatory integrals. IMA Journal of Numerical Analysis 31 (2011), no 4, 1253-1280. Preprint available http://opus.bath.ac.uk/26655

The algorithms for the implementation are detailed and fully explained there.

In short, we deal with the approximation of $\displaystyle\int_{a}^b f(t)\exp({\rm i}k t)\,{\rm d}t$ $\displaystyle\int_{-1}^{1} f(t)\log((t-\alpha)^2)\exp({\rm i}k t)\,{\rm d}t$

(Notice that in the second integral the interval is fixed). Here $k$ is real and $\alpha\in[-1,1]$. Actually $k=0$ and we have the Classical Clenshaw-Curtis rule in the first case and a product Clenshaw-Curtis rule in the second case.

These rules are based on replacing $f$ by the interpolating polynomial at the Chebyshev nodes. Hence, in the first case, and with $a=-1,\ b=1$ $\displaystyle\int_{-1}^1 f(t)\exp({\rm i}k t)\,{\rm d}t \approx \int_{-1}^1 p_n(t)\exp({\rm i}k t)\,{\rm d}t$

where $\displaystyle \mathbb{P}_n \ni p_n \quad \text{s.t.}\quad p_n\left(\cos\frac{\pi j}{n} \right)=f\left(\cos\frac{\pi j}{n} \right),\quad j=0,\ldots,n$

The second integral is computed exactly. This can be done, in theory. In practise the analytical expression are both complicate and unstable as $n$ increases. The implementation is based on the following lines:

1. Compute  $latex p_n$ in the Chebyshev basis $\{T_m\}_{m=0}^n$ , which can be done efficiently using FFT techniques
2. Compute $\displaystyle \int_{-1}^1 T_m(s)\exp(iks)\,{\rm d}s$

via a fast recursive relations which are shown to be stable.

The other integral is dealt following similar ideas.

You can download an updated version, as of July 2017, of the Matlab implementation of such algorithms. Fejer-type quadrature rules are supported as well (not analysed in the cited papers, though)

I appreciate any comment, suggestion or (even much more) correction you can share with me.