Fractional Laplacian for generalized Zernike Disk polynomials #73

TSGut · 2021-05-03T10:32:01Z

Starting work on fractional Laplacian here. Should also include an example file for fractional Poisson by the time this is done.

TSGut · 2021-05-03T10:34:42Z

Might be a good time to discuss what we want the syntax to look like for this object @dlfivefifty. At the moment I just have FractionalLaplacian(β) to represent (-Δ)^(β) acting on (1-r^2)^β * Zernike(β).

codecov · 2021-05-03T10:39:24Z

Codecov Report

Merging #73 (4a1b886) into master (097dc70) will increase coverage by 0.07%.
The diff coverage is 95.00%.

@@            Coverage Diff             @@
##           master      #73      +/-   ##
==========================================
+ Coverage   94.73%   94.81%   +0.07%     
==========================================
  Files           2        3       +1     
  Lines         456      482      +26     
==========================================
+ Hits          432      457      +25     
- Misses         24       25       +1

Impacted Files	Coverage Δ
src/disk.jl	`98.70% <95.00%> (-0.67%)`	⬇️
src/ModalInterlace.jl	`100.00% <0.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 097dc70...4a1b886. Read the comment docs.

dlfivefifty · 2021-05-10T09:08:50Z

I would say we support either abs(Δ)^α or (-Δ)^α where Δ = Laplacian(axes(Z,1)).

Maybe add AbsLaplacianPower:

struct AbsLaplacianPower{T,D,A} <: LazyQuasiMatrix{T}
    axis::Inclusion{T,D}
    α::A
end

AbsLaplacianPower{T}(axis::Inclusion{<:Any,D},α) where {T,D} = Laplacian{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(Δ,1)
-(Δ::Laplacian) = abs(Δ)

^(D::AbsLaplacianPower, k) = AbsLaplacianPower(D.axis, D.α*k)

Probably add it here https://github.com/JuliaApproximation/HarmonicOrthogonalPolynomials.jl/blob/master/src/laplace.jl

merge updates

TSGut · 2021-05-15T11:54:48Z

@dlfivefifty Okay, I think this is basically doing everything I wanted now. Fractional Laplacians and fractional Poisson equation on disk work and pass various explicit solution tests. This is based on results in https://link.springer.com/article/10.1007%2Fs00365-016-9336-4 and similar ideas seem to have been explored before for spectral methods, see e.g. https://arxiv.org/abs/1812.08325. We can now adaptively do all of the stuff in that paper in 2D and it would be easy to adapt for arbitrary dimensions if we had the OPs set up for that.

I was thinking we could do nice surface plots for the fractional Poisson equation but plotting Zernike and weighted Zernike via surface() seems broken? Is disk plotting functional yet and do we have a minimal working example somewhere?

dlfivefifty · 2021-05-20T18:37:54Z

Plotting should work with plotly() backend

TSGut · 2021-05-22T13:51:28Z

Yes but it seems I can only get it to plot things in Zernike basis, but get errors for weighted bases. Are weighted bases also meant to be working?

julia> using MultivariateOrthogonalPolynomials, InfiniteArrays, Plots

julia> plotly()

julia> Z = Zernike(1)
Zernike{Float64}

julia> W = Weighted(Z)
Weighted(Zernike{Float64})

julia> u = Z*Vcat(1,zeros(∞))
Zernike{Float64} * [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  ]

julia> w = W*Vcat(1,zeros(∞))
Weighted(Zernike{Float64}) * [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  ]

julia> surface(u) # this works as expected

julia> surface(w)
ERROR: LoadError: Overload Grid
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] _grid(#unused#::ContinuumArrays.BasisLayout, P::Zernike{Float64})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:124
  [3] grid(P::Zernike{Float64})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:129
  [4] _grid(#unused#::ContinuumArrays.SubBasisLayout, P::QuasiArrays.SubQuasiArray{Float64, 2, Zernike{Float64}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:127
  [5] grid(P::QuasiArrays.SubQuasiArray{Float64, 2, Zernike{Float64}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:129
  [6] _grid(#unused#::ContinuumArrays.SubWeightedBasisLayout, P::QuasiArrays.SubQuasiArray{Float64, 2, Weighted{Float64, Zernike{Float64}}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:128
  [7] grid(P::QuasiArrays.SubQuasiArray{Float64, 2, Weighted{Float64, Zernike{Float64}}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/bases/bases.jl:129
  [8] _plotgrid(#unused#::ContinuumArrays.SubWeightedBasisLayout, P::QuasiArrays.SubQuasiArray{Float64, 2, Weighted{Float64, Zernike{Float64}}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:11
  [9] plotgrid(g::QuasiArrays.SubQuasiArray{Float64, 2, Weighted{Float64, Zernike{Float64}}, Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitBall{StaticArrays.SVector{2, Float64}, :closed}}, UnitRange{Int64}}, false})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:13
 [10] _mul_plotgrid(#unused#::Tuple{ContinuumArrays.WeightedBasisLayout, LazyArrays.PaddedLayout{LazyArrays.ApplyLayout{typeof(vcat)}}}, ::Tuple{Weighted{Float64, Zernike{Float64}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Int64, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:4
 [11] _plotgrid(lay::LazyArrays.ApplyLayout{typeof(*)}, P::QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Weighted{Float64, Zernike{Float64}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Int64, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:8
 [12] plotgrid(g::QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Weighted{Float64, Zernike{Float64}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Int64, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:13
 [13] macro expansion
    @ ~/.julia/packages/ContinuumArrays/6hwP3/src/plotting.jl:19 [inlined]
 [14] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, g::QuasiArrays.AbstractQuasiVector{T} where T)
    @ ContinuumArrays ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:282
 [15] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/CirY4/src/user_recipe.jl:36
 [16] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/CirY4/src/RecipesPipeline.jl:70
 [17] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/vVVub/src/plot.jl:172
 [18] plot(args::Any; kw::Any)
    @ Plots ~/.julia/packages/Plots/vVVub/src/plot.jl:58
 [19] surface(args::Any; kw::Any)
    @ Plots ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:403
 [20] surface(args::Any)
    @ Plots ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:403
 [21] top-level scope
    @ ~/Documents/Projects/MultivariateOrthogonalPolynomials.jl/examples/diskfractionalpoisson.jl:92

dlfivefifty · 2021-05-23T12:36:55Z

Ah yes that's a missing feature...probably need to add to ContinuumArrays.jl/src/plotting.jl something like

_plotgrid(::WeightedBasisLayouts, wP) = plotgrid(unweightedbasis(wP))

Do you feel like trying to make that change and see if it works?

(We do need to add some more features to get plotting to take advantage of fast transforms but that I can wait.)

TSGut · 2021-05-23T17:10:02Z

Yes, I'll have a look! 👍 I'm curious to see a bit of the plotting backend, haven't had too much of a reason to look into it before.

dlfivefifty · 2021-05-23T17:20:53Z

It’s very simple like 10 lines

TSGut · 2021-05-24T08:41:04Z

A lot simpler than I thought it would be. And yes it looks like that works. I'll do some more tests and then add a PR to ContinuumArrays.jl.

src/disk.jl

TSGut · 2021-05-24T12:22:40Z

With JuliaApproximation/ContinuumArrays.jl#96 we can now also do some nice plots for the fractional Poisson equation on the disk.

TSGut added 2 commits May 3, 2021 11:28

Start work on fractional Laplacian on generalized Zernike

248a629

standard Laplacian special case test

e07dd0a

TSGut added 7 commits May 3, 2021 11:52

change convention and assert basis

d9fb00e

prep for merging of updates

c899751

prep

c15d7dc

update files

4349a81

Merge branch 'master' into master

8c56090

tests

5c6a26b

test fix imports

3664cfe

TSGut changed the title ~~WIP: Fractional Laplacian for generalized Zernike DIsk polynomials~~ WIP: Fractional Laplacian for generalized Zernike Disk polynomials May 10, 2021

TSGut mentioned this pull request May 10, 2021

Powers of negative Laplacians JuliaApproximation/HarmonicOrthogonalPolynomials.jl#32

Merged

TSGut added 4 commits May 13, 2021 14:47

Merge pull request #2 from JuliaApproximation/master

8db1743

merge updates

update to AbsLaplacianPower

a201c75

tests for fractional laplacian: constant, radial and non radial

bb6a04e

Tests for explicit solutions of Fractional Poisson Eq.

3346636

TSGut mentioned this pull request May 24, 2021

Plotting for weighted bases JuliaApproximation/ContinuumArrays.jl#96

Merged

dlfivefifty approved these changes May 24, 2021

View reviewed changes

src/disk.jl Outdated Show resolved Hide resolved

TSGut added 2 commits May 24, 2021 13:01

fix misleading comments

6533f26

tweak for BigFloat support

4a1b886

TSGut changed the title ~~WIP: Fractional Laplacian for generalized Zernike Disk polynomials~~ Fractional Laplacian for generalized Zernike Disk polynomials May 24, 2021

dlfivefifty merged commit 04704a8 into JuliaApproximation:master May 24, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Fractional Laplacian for generalized Zernike Disk polynomials #73

Fractional Laplacian for generalized Zernike Disk polynomials #73

Uh oh!

TSGut commented May 3, 2021

Uh oh!

TSGut commented May 3, 2021 •

edited

Loading

Uh oh!

codecov bot commented May 3, 2021 •

edited

Loading

Uh oh!

dlfivefifty commented May 10, 2021

Uh oh!

TSGut commented May 15, 2021 •

edited

Loading

Uh oh!

dlfivefifty commented May 20, 2021

Uh oh!

TSGut commented May 22, 2021 •

edited

Loading

Uh oh!

dlfivefifty commented May 23, 2021 •

edited

Loading

Uh oh!

TSGut commented May 23, 2021

Uh oh!

dlfivefifty commented May 23, 2021

Uh oh!

TSGut commented May 24, 2021

Uh oh!

Uh oh!

TSGut commented May 24, 2021

Uh oh!

Uh oh!

Fractional Laplacian for generalized Zernike Disk polynomials #73

Fractional Laplacian for generalized Zernike Disk polynomials #73

Uh oh!

Conversation

TSGut commented May 3, 2021

Uh oh!

TSGut commented May 3, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov bot commented May 3, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

dlfivefifty commented May 10, 2021

Uh oh!

TSGut commented May 15, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dlfivefifty commented May 20, 2021

Uh oh!

TSGut commented May 22, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dlfivefifty commented May 23, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

TSGut commented May 23, 2021

Uh oh!

dlfivefifty commented May 23, 2021

Uh oh!

TSGut commented May 24, 2021

Uh oh!

Uh oh!

TSGut commented May 24, 2021

Uh oh!

Uh oh!

TSGut commented May 3, 2021 •

edited

Loading

codecov bot commented May 3, 2021 •

edited

Loading

TSGut commented May 15, 2021 •

edited

Loading

TSGut commented May 22, 2021 •

edited

Loading

dlfivefifty commented May 23, 2021 •

edited

Loading