Нелокальная диффузия в S^2

В этом примере вычисляется спектр оператора нелокальной диффузии:

определенного в Eq. (2) в работе

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.

В приведенном выше примере , и ядро:

где — индикаторная функция в множестве .

Этот нелокальный оператор диагонализируется сферическими гармониками:

и его собственные функции задаются обобщенной формулой Функа-Гекке:

В статье авторы используют квадратуру Кленшо-Кертиса и асимптотическую оценку многочленов Лежандра для достижения сложности для оценки первых собственных значений. При смене базиса эта сложность может быть снижена до .

Сначала представим:

Затем представим с многочленами Якоби и проинтегрируем с помощью DLMF 18.9.16:

Приведенный ниже код реализует этот алгоритм, используя преобразование Якоби-Якоби plan_jac2jac. Для обеспечения численной стабильности преобразование многочленов Якоби в делится на преобразование из в , а затем на преобразование из в .

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)

Спектр в 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

Спектр в 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

Относительная погрешность -нормы:

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

0.0