Skip to content

Commit 04704a8

Browse files
TSGutTSGut
authored
Fractional Laplacian for generalized Zernike Disk polynomials (#73)
* Start work on fractional Laplacian on generalized Zernike * standard Laplacian special case test * change convention and assert basis * prep for merging of updates * prep * update files * tests * test fix imports * update to AbsLaplacianPower * tests for fractional laplacian: constant, radial and non radial * Tests for explicit solutions of Fractional Poisson Eq. * fix misleading comments * tweak for BigFloat support
1 parent 793dd20 commit 04704a8

File tree

4 files changed

+285
-5
lines changed

4 files changed

+285
-5
lines changed

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, Multivariat
2424

2525
export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian,
2626
MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
27-
zerniker, zernikez, Weighted, Block, ZernikeWeight
27+
zerniker, zernikez, Weighted, Block, ZernikeWeight, AbsLaplacianPower
2828

2929
include("ModalInterlace.jl")
3030
include("disk.jl")

‎src/disk.jl

Lines changed: 40 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -237,7 +237,7 @@ function Base.view(W::WeightedZernikeLaplacianDiag{T}, K::Block{1}) where T
237237
convert(AbstractVector{T}, -4*interlace(v, v[2:end]))
238238
else
239239
v = (d:K̃) .* (d-1:(-1):1)
240-
convert(AbstractVector{T}, -4 * interlace(v, v))
240+
convert(AbstractVector{T}, -4*interlace(v, v))
241241
end
242242
end
243243

@@ -248,6 +248,45 @@ getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,
248248
WZ.P * Diagonal(WeightedZernikeLaplacianDiag{eltype(eltype(WZ))}())
249249
end
250250

251+
###
252+
# Fractional Laplacian
253+
###
254+
255+
function *(L::AbsLaplacianPower, WZ::Weighted{<:Any,<:Zernike{<:Any}})
256+
@assert axes(L,1) == axes(WZ,1) && WZ.P.a == 0 && WZ.P.b == L.α
257+
WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(L.α)}(L.α))
258+
end
259+
260+
# gives the entries for the (negative!) fractional Laplacian (-Δ)^(α) times (1-r^2)^α * Zernike(α)
261+
struct WeightedZernikeFractionalLaplacianDiag{T} <: AbstractBlockVector{T}
262+
α::T
263+
end
264+
265+
axes(::WeightedZernikeFractionalLaplacianDiag) = (blockedrange(oneto(∞)),)
266+
copy(R::WeightedZernikeFractionalLaplacianDiag) = R
267+
268+
MemoryLayout(::Type{<:WeightedZernikeFractionalLaplacianDiag}) = LazyLayout()
269+
Base.BroadcastStyle(::Type{<:Diagonal{<:Any,<:WeightedZernikeFractionalLaplacianDiag}}) = LazyArrayStyle{2}()
270+
271+
getindex(W::WeightedZernikeFractionalLaplacianDiag, k::Integer) = W[findblockindex(axes(W,1),k)]
272+
273+
function Base.view(W::WeightedZernikeFractionalLaplacianDiag{T}, K::Block{1}) where T
274+
l = Int(K)
275+
if isodd(l)
276+
m = Vcat(0,interlace(Array(2:2:l),Array(2:2:l)))
277+
else #if iseven(l)
278+
m = Vcat(interlace(Array(1:2:l),Array(1:2:l)))
279+
end
280+
return convert(AbstractVector{T}, 2^(2*W.α)*fractionalcfs2d.(l-1,m,W.α))
281+
end
282+
283+
# generic d-dimensional ball fractional coefficients without the 2^(2*β) factor. m is assumed to be entered as abs(m)
284+
function fractionalcfs(l::Integer, m::Integer, α::T, d::Integer) where T
285+
n = (l-m)÷2
286+
return exp(loggamma(α+n+1)+loggamma((2*α+2*n+d+2*m)/2)-loggamma(one(T)+n)-loggamma((2*one(T)*m+2*n+d)/2))
287+
end
288+
# 2 dimensional special case, again without the 2^(2*β) factor
289+
fractionalcfs2d(l::Integer, m::Integer, β) = fractionalcfs(l,m,β,2)
251290

252291
function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
253292
TV = promote_type(T,V)

‎test/test_disk.jl

Lines changed: 243 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test
2-
import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform
1+
using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test, SpecialFunctions
2+
import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeTransform, ZernikeITransform, *
33
import ClassicalOrthogonalPolynomials: HalfWeighted
44

55

@@ -245,3 +245,244 @@ import ClassicalOrthogonalPolynomials: HalfWeighted
245245
@test rep[1].args == (first.(g),last.(g),u[g])
246246
end
247247
end
248+
249+
@testset "Fractional Laplacian on Unit Disk" begin
250+
@testset "Fractional Laplacian on Disk: (-Δ)^(β) == -Δ when β=1" begin
251+
WZ = Weighted(Zernike(1.))
252+
Δ = Laplacian(axes(WZ,1))
253+
Δ_Z = Zernike(1) \ (Δ * WZ)
254+
Δfrac = AbsLaplacianPower(axes(WZ,1),1.)
255+
Δ_Zfrac = Zernike(1) \ (Δfrac * WZ)
256+
@test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100]
257+
end
258+
259+
@testset "Fractional Laplacian on Disk: Computing f where (-Δ)^(β) u = f" begin
260+
@testset "Set 1 - Explicitly known constant f" begin
261+
# set up basis
262+
β = 1.34
263+
Z = Zernike(β)
264+
WZ = Weighted(Z)
265+
xy = axes(WZ,1)
266+
x,y = first.(xy),last.(xy)
267+
# generate fractional Laplacian
268+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
269+
Δ_Zfrac = Z \ (Δfrac * WZ)
270+
# define function whose fractional Laplacian is known
271+
u = @. (1 - x^2 - y^2).^β
272+
# explicit and computed solutions
273+
fexplicit0(d,α) = 2^α*gamma(α/2+1)*gamma((d+α)/2)/gamma(d/2) # note that here, α = 2*β
274+
f = Z*(Δ_Zfrac*(WZ \ u))
275+
# compare
276+
@test fexplicit0(2,2*β) ≈ f[(0.1,0.4)] ≈ f[(0.1137,0.001893)] ≈ f[(0.3721,0.3333)]
277+
278+
# again for different β
279+
β = 2.11
280+
Z = Zernike(β)
281+
WZ = Weighted(Z)
282+
xy = axes(WZ,1)
283+
x,y = first.(xy),last.(xy)
284+
# generate fractional Laplacian
285+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
286+
Δ_Zfrac = Z \ (Δfrac * WZ)
287+
# define function whose fractional Laplacian is known
288+
u = @. (1 - x^2 - y^2).^β
289+
# computed solution
290+
f = Z*(Δ_Zfrac*(WZ \ u))
291+
# compare
292+
@test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1731,0.091893)] ≈ f[(0.3791,0.333333)]
293+
294+
# again for different β
295+
β = 3.14159
296+
Z = Zernike(β)
297+
WZ = Weighted(Z)
298+
xy = axes(WZ,1)
299+
x,y = first.(xy),last.(xy)
300+
# generate fractional Laplacian
301+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
302+
Δ_Zfrac = Z \ (Δfrac * WZ)
303+
# define function whose fractional Laplacian is known
304+
u = @. (1 - x^2 - y^2).^β
305+
# computed solution
306+
f = Z*(Δ_Zfrac*(WZ \ u))
307+
# compare
308+
@test fexplicit0(2,2*β) ≈ f[(0.14,0.41)] ≈ f[(0.1837,0.101893)] ≈ f[(0.37222,0.2222)]
309+
end
310+
@testset "Set 2 - Explicitly known radially symmetric f" begin
311+
β = 1.1
312+
Z = Zernike(β)
313+
WZ = Weighted(Z)
314+
xy = axes(WZ,1)
315+
x,y = first.(xy),last.(xy)
316+
# generate fractional Laplacian
317+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
318+
Δ_Zfrac = Z \ (Δfrac * WZ)
319+
# define function whose fractional Laplacian is known
320+
u = @. (1 - x^2 - y^2).^(β+1)
321+
# explicit and computed solutions
322+
fexplicit1(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
323+
f = Z*(Δ_Zfrac*(WZ \ u))
324+
# compare
325+
@test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
326+
@test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
327+
@test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
328+
329+
# again for different β
330+
β = 2.71999
331+
Z = Zernike(β)
332+
WZ = Weighted(Z)
333+
xy = axes(WZ,1)
334+
x,y = first.(xy),last.(xy)
335+
# generate fractional Laplacian
336+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
337+
Δ_Zfrac = Z \ (Δfrac * WZ)
338+
# define function whose fractional Laplacian is known
339+
u = @. (1 - x^2 - y^2).^(β+1)
340+
# explicit and computed solutions
341+
f = Z*(Δ_Zfrac*(WZ \ u))
342+
# compare
343+
@test fexplicit1(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
344+
@test fexplicit1(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
345+
@test fexplicit1(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
346+
end
347+
@testset "Set 3 - Explicitly known f, not radially symmetric" begin
348+
# dependence on x
349+
β = 2.71999
350+
Z = Zernike(β)
351+
WZ = Weighted(Z)
352+
xy = axes(WZ,1)
353+
x,y = first.(xy),last.(xy)
354+
# generate fractional Laplacian
355+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
356+
Δ_Zfrac = Z \ (Δfrac * WZ)
357+
# define function whose fractional Laplacian is known
358+
u = @. (1 - x^2 - y^2).^(β)*x
359+
# explicit and computed solutions
360+
fexplicit2(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[1] # α = 2*β
361+
f = Z*(Δ_Zfrac*(WZ \ u))
362+
# compare
363+
@test fexplicit2(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
364+
@test fexplicit2(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
365+
@test fexplicit2(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
366+
367+
# different β, dependence on y
368+
β = 1.91239
369+
Z = Zernike(β)
370+
WZ = Weighted(Z)
371+
xy = axes(WZ,1)
372+
x,y = first.(xy),last.(xy)
373+
# generate fractional Laplacian
374+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
375+
Δ_Zfrac = Z \ (Δfrac * WZ)
376+
# define function whose fractional Laplacian is known
377+
u = @. (1 - x^2 - y^2).^(β)*y
378+
# explicit and computed solutions
379+
fexplicit3(d,α,x) = 2^α*gamma(α/2+1)*gamma((d+α)/2+1)/gamma(d/2+1)*x[2] # α = 2*β
380+
f = Z*(Δ_Zfrac*(WZ \ u))
381+
# compare
382+
@test fexplicit3(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
383+
@test fexplicit3(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
384+
@test fexplicit3(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
385+
end
386+
@testset "Set 4 - Explicitly known f, different non-radially-symmetric example" begin
387+
# dependence on x
388+
β = 1.21999
389+
Z = Zernike(β)
390+
WZ = Weighted(Z)
391+
xy = axes(WZ,1)
392+
x,y = first.(xy),last.(xy)
393+
# generate fractional Laplacian
394+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
395+
Δ_Zfrac = Z \ (Δfrac * WZ)
396+
# define function whose fractional Laplacian is known
397+
u = @. (1 - x^2 - y^2).^(β+1)*x
398+
# explicit and computed solutions
399+
fexplicit4(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
400+
f = Z*(Δ_Zfrac*(WZ \ u))
401+
# compare
402+
@test fexplicit4(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
403+
@test fexplicit4(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
404+
@test fexplicit4(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
405+
406+
# different β, dependence on y
407+
β = 0.141
408+
Z = Zernike(β)
409+
WZ = Weighted(Z)
410+
xy = axes(WZ,1)
411+
x,y = first.(xy),last.(xy)
412+
# generate fractional Laplacian
413+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
414+
Δ_Zfrac = Z \ (Δfrac * WZ)
415+
# define function whose fractional Laplacian is known
416+
u = @. (1 - x^2 - y^2).^(β+1)*y
417+
# explicit and computed solutions
418+
fexplicit5(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
419+
f = Z*(Δ_Zfrac*(WZ \ u))
420+
# compare
421+
@test fexplicit5(2,2*β,(0.94,0.01)) ≈ f[(0.94,0.01)]
422+
@test fexplicit5(2,2*β,(0.14,0.41)) ≈ f[(0.14,0.41)]
423+
@test fexplicit5(2,2*β,(0.221,0.333)) ≈ f[(0.221,0.333)]
424+
end
425+
426+
@testset "Fractional Poisson equation on Disk: Comparison with explicitly known solutions" begin
427+
@testset "Set 1 - Radially symmetric solution" begin
428+
# define basis
429+
β = 1.1812
430+
Z = Zernike(β)
431+
WZ = Weighted(Z)
432+
xy = axes(WZ,1)
433+
x,y = first.(xy),last.(xy)
434+
# generate fractional Laplacian
435+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
436+
Δ_Zfrac = Z \ (Δfrac * WZ)
437+
# define function whose fractional Laplacian is known
438+
uexplicit = @. (1 - x^2 - y^2).^(β+1)
439+
uexplicitcfs = WZ \ uexplicit
440+
# RHS
441+
RHS(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2)/gamma(d/2)*(1-(1+α/d)*norm(x)^2) # α = 2*β
442+
RHScfs = Z \ @. RHS.(2,2*β,xy)
443+
# compute solution
444+
ucomputed = Δ_Zfrac \ RHScfs
445+
@test uexplicitcfs[1:100] ≈ ucomputed[1:100]
446+
end
447+
@testset "Set 2 - Non-radially-symmetric solutions" begin
448+
# dependence on y
449+
β = 0.98812
450+
Z = Zernike(β)
451+
WZ = Weighted(Z)
452+
xy = axes(WZ,1)
453+
x,y = first.(xy),last.(xy)
454+
# generate fractional Laplacian
455+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
456+
Δ_Zfrac = Z \ (Δfrac * WZ)
457+
# define function whose fractional Laplacian is known
458+
uexplicit = @. (1 - x^2 - y^2).^(β+1)*y
459+
uexplicitcfs = WZ \ uexplicit
460+
# RHS
461+
RHS2(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[2] # α = 2*β
462+
RHS2cfs = Z \ @. RHS2.(2,2*β,xy)
463+
# compute solution
464+
ucomputed = Δ_Zfrac \ RHS2cfs
465+
@test uexplicitcfs[1:100] ≈ ucomputed[1:100]
466+
467+
# different β, dependence on x
468+
β = 0.506
469+
Z = Zernike(β)
470+
WZ = Weighted(Z)
471+
xy = axes(WZ,1)
472+
x,y = first.(xy),last.(xy)
473+
# generate fractional Laplacian
474+
Δfrac = AbsLaplacianPower(axes(WZ,1),β)
475+
Δ_Zfrac = Z \ (Δfrac * WZ)
476+
# define function whose fractional Laplacian is known
477+
uexplicit = @. (1 - x^2 - y^2).^(β+1)*x
478+
uexplicitcfs = WZ \ uexplicit
479+
# RHS
480+
RHS3(d,α,x) = 2^α*gamma(α/2+2)*gamma((d+α)/2+1)/gamma(d/2+1)*(1-(1+α/(d+2))*norm(x)^2)*x[1] # α = 2*β
481+
RHS3cfs = Z \ @. RHS3.(2,2*β,xy)
482+
# compute solution
483+
ucomputed = Δ_Zfrac \ RHS3cfs
484+
@test uexplicitcfs[1:100] ≈ ucomputed[1:100]
485+
end
486+
end
487+
end
488+
end

‎test/test_domaindecomposition.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ using ApproxFun, MultivariateOrthogonalPolynomials, Plots, BlockArrays
55
# Neumann
66

77
d = Triangle()
8-
S = DirichletTriangle{1,1,1}(d))
8+
S = DirichletTriangle{1,1,1}(d)
99
Dx = Derivative(S,[1,0])
1010
Dy = Derivative(S,[0,1])
1111
B₁ = I : S → Legendre(Vec(0,0) .. Vec(0,1))

0 commit comments

Comments
 (0)