Numerical computation of the half Laplacian by means of a fast convolution algorithm
share
In this paper, we develop a fast and accurate pseudospectral method to approximate numerically the fractional Laplacian (-Δ)^α/2 of a function on ℝ for α=1; this case, commonly referred to as the half Laplacian, is equivalent to the Hilbert transform of the derivative of the function. The main ideas are as follows. Given a twice continuously differentiable bounded function u∈𝐶_b^2(ℝ), we apply the change of variable x=L(s), with L>0 and s∈[0,π], which maps ℝ into [0,π], and denote (-Δ)_s^1/2u(x(s)) ≡ (-Δ)^1/2u(x). Therefore, by performing a Fourier series expansion of u(x(s)), the problem is reduced to computing (-Δ)_s^1/2e^iks≡ (-Δ)^1/2(x + i)^k/(1+x^2)^k/2. On a previous work, we considered the case with k even and α∈(0,2), so we focus now on the case with k odd. More precisely, we express (-Δ)_s^1/2e^iks for k odd in terms of the Gaussian hypergeometric function _2F_1, and also as a well-conditioned finite sum. Then, we use a fast convolution result, that enable us to compute very efficiently ∑_l = 0^Ma_l(-Δ)_s^1/2e^i(2l+1)s, for extremely large values of M. This enables us to approximate (-Δ)_s^1/2u(x(s)) in a fast and accurate way, especially when u(x(s)) is not periodic of period π.
READ FULL TEXT