Powers of negative Laplacians (#32)

* import export

* abs power laplacian struct

* tests

* v0.1

Author: TSGut

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicOrthogonalPolynomials"
 uuid = "e416a80e-9640-42f3-8df8-80a93ca01ea5"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.0.5"
+version = "0.1"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"

src/HarmonicOrthogonalPolynomials.jl

@@ -1,7 +1,7 @@
 module HarmonicOrthogonalPolynomials
 using FastTransforms, LinearAlgebra, ClassicalOrthogonalPolynomials, ContinuumArrays, DomainSets, 
         BlockArrays, BlockBandedMatrices, InfiniteArrays, StaticArrays, QuasiArrays, Base, SpecialFunctions
-import Base: OneTo, axes, getindex, convert, to_indices, _maybetail, tail, eltype, *, ==, ^, copy
+import Base: OneTo, axes, getindex, convert, to_indices, _maybetail, tail, eltype, *, ==, ^, copy, -, abs
 import BlockArrays: block, blockindex, unblock, BlockSlice
 import DomainSets: indomain
 import LinearAlgebra: norm, factorize
@@ -12,7 +12,7 @@ import BlockBandedMatrices: BlockRange1
 import FastTransforms: Plan, interlace
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
 
-export SphericalHarmonic, UnitSphere, SphericalCoordinate, RadialCoordinate, Block, associatedlegendre, RealSphericalHarmonic, sphericalharmonicy, Laplacian
+export SphericalHarmonic, UnitSphere, SphericalCoordinate, RadialCoordinate, Block, associatedlegendre, RealSphericalHarmonic, sphericalharmonicy, Laplacian, AbsLaplacianPower, abs, -, ^
 
 include("multivariateops.jl")
 include("spheretrav.jl")

src/laplace.jl

@@ -1,3 +1,5 @@
+# Laplacian
+
 struct Laplacian{T,D} <: LazyQuasiMatrix{T}
     axis::Inclusion{T,D}
 end
@@ -13,4 +15,22 @@ copy(D::Laplacian) = Laplacian(copy(D.axis))
 @simplify function *(Δ::Laplacian, P::AbstractSphericalHarmonic)
      # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
      P * Diagonal(mortar(Fill.((-(0:∞)-(0:∞).^2), 1:2:∞)))
-end
+end
+
+# Negative fractional Laplacian (-Δ)^α or equiv. abs(Δ)^α
+
+struct AbsLaplacianPower{T,D,A} <: LazyQuasiMatrix{T}
+    axis::Inclusion{T,D}
+    α::A
+end
+
+AbsLaplacianPower{T}(axis::Inclusion{<:Any,D},α) where {T,D} = AbsLaplacianPower{T,D,typeof(α)}(axis,α)
+
+axes(D:: AbsLaplacianPower) = (D.axis, D.axis)
+==(a:: AbsLaplacianPower, b:: AbsLaplacianPower) = a.axis == b.axis && a.α == b.α
+copy(D:: AbsLaplacianPower) = AbsLaplacianPower(copy(D.axis), D.α)
+
+abs(Δ::Laplacian) = AbsLaplacianPower(axes(Δ,1),1)
+-(Δ::Laplacian) = abs(Δ)
+
+^(D::AbsLaplacianPower, k) = AbsLaplacianPower(D.axis, D.α*k)

test/runtests.jl

@@ -374,3 +374,30 @@ end
     @test (Δ^2*S*(S\f1.(xyz)))[t] ≈ (Δ*Δ*S*(S\f1.(xyz)))[t] ≈ Δ2_f1(t)
     @test (Δ^3*S*(S\f1.(xyz)))[t] ≈ (Δ*Δ*Δ*S*(S\f1.(xyz)))[t] ≈ Δ3_f1(t)
 end
+
+@testset "abs(Δ)^α - Basics of absolute Laplacian powers" begin
+    # Set 1
+    α = 1/3
+    S = SphericalHarmonic()
+    Sxyz = axes(S,1)
+    SΔα = AbsLaplacianPower(Sxyz,α)
+    Δ = Laplacian(Sxyz)
+    @test copy(SΔα) == SΔα
+    @test SΔα isa AbsLaplacianPower
+    @test SΔα isa QuasiArrays.LazyQuasiMatrix
+    @test axes(SΔα) == (axes(S,1),axes(S,1))
+    @test abs(Δ) == -Δ == AbsLaplacianPower(axes(Δ,1),1)
+    @test abs(Δ)^α == SΔα
+    # Set 2
+    α = 7/13
+    S = SphericalHarmonic()
+    Sxyz = axes(S,1)
+    SΔα = AbsLaplacianPower(Sxyz,α)
+    Δ = Laplacian(Sxyz)
+    @test copy(SΔα) == SΔα
+    @test SΔα isa AbsLaplacianPower
+    @test SΔα isa QuasiArrays.LazyQuasiMatrix
+    @test axes(SΔα) == (axes(S,1),axes(S,1))
+    @test abs(Δ) == -Δ == AbsLaplacianPower(axes(Δ,1),1)
+    @test abs(Δ)^α == SΔα
+end