Nonlocal diffusion on S^2

This example calculates the spectrum of the nonlocal diffusion operator:

defined in Eq. (2) of

R.M. Slevinsky, H. Montanelli, and Q. Du, A spectral method for nonlocal diffusion operators on the sphere, J. Comp. Phys., 372:893—911, 2018.

In the above, , , and the kernel:

where is the indicator function on the set .

This nonlocal operator is diagonalized by spherical harmonics:

and its eigenfunctions are given by the generalized Funk—Hecke formula:

In the paper, the authors use Clenshaw—Curtis quadrature and asymptotic evaluation of Legendre polynomials to achieve complexity for the evaluation of the first eigenvalues. With a change of basis, this complexity can be reduced to .

First, we represent:

Then, we represent with Jacobi polynomials and we integrate using DLMF 18.9.16:

The code below implements this algorithm, making use of the Jacobi—Jacobi transform plan_jac2jac. For numerical stability, the conversion from Jacobi polynomials to is divided into conversion from to , before conversion from to .

using FastTransforms, LinearAlgebra

function oprec!(n::Integer, v::AbstractVector, alpha::Real, delta2::Real)
    if n > 0
        v[1] = 1
    end
    if n > 1
        v[2] = (4*alpha+8-(alpha+4)*delta2)/4
    end
    for i = 1:n-2
        v[i+2] = (((2*i+alpha+2)*(2*i+alpha+4)+alpha*(alpha+2))/(2*(i+1)*(2*i+alpha+2))*(2*i+alpha+3)/(i+alpha+3) - delta2/4*(2*i+alpha+3)/(i+1)*(2*i+alpha+4)/(i+alpha+3))*v[i+1] - (i+alpha+1)/(i+alpha+3)*(2*i+alpha+4)/(2*i+alpha+2)*v[i]
    end
    return v
end

function evaluate_lambda(n::Integer, alpha::T, delta::T) where T
    delta2 = delta*delta
    scl = (1+alpha)*(2-delta2/2)

    lambda = Vector{T}(undef, n)

    if n > 0
        lambda[1] = 0
    end
    if n > 1
        lambda[2] = -2
    end

    oprec!(n-2, view(lambda, 3:n), alpha, delta2)

    for i = 2:n-1
        lambda[i+1] *= -scl/(i-1)
    end

    p = plan_jac2jac(T, n-1, zero(T), zero(T), alpha, zero(T))

    lmul!(p', view(lambda, 2:n))

    for i = 2:n-1
        lambda[i+1] = ((2i-1)*lambda[i+1] + (i-1)*lambda[i])/i
    end

    for i = 2:n-1
        lambda[i+1] += lambda[i]
    end

    return lambda
end

evaluate_lambda (generic function with 1 method)

The spectrum in Float64:

lambda = evaluate_lambda(10, -0.5, 1.0)

10-element Vector{Float64}:
   0.0
  -2.0
  -5.5
  -9.75
 -14.09375
 -18.203125
 -22.08984375
 -25.935546875
 -29.8870849609375
 -33.95416259765625

The spectrum in BigFloat:

lambdabf = evaluate_lambda(10, parse(BigFloat, "-0.5"), parse(BigFloat, "1.0"))

10-element Vector{BigFloat}:
   0.0
  -2.0
  -5.5
  -9.75
 -14.09375
 -18.203125
 -22.08984375
 -25.935546875
 -29.8870849609375
 -33.95416259765625

The -norm relative error:

norm(lambda-lambdabf, Inf)/norm(lambda, Inf)

0.0

This page was generated using Literate.jl.