Example of using p-FEM (#173)

* Example of using p-FEM

* Fix bugs

* update tests

* Update index.md

* v0.12.6

Author: dlfivefifty

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.12.5"
+version = "0.12.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "1.3.1"
 BandedMatrices = "0.17.17, 1"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.12"
-ContinuumArrays = "0.17"
+ContinuumArrays = "0.17.3"
 DomainSets = "0.6, 0.7"
 FFTW = "1.1"
 FastGaussQuadrature = "0.5, 1"

docs/src/index.md

@@ -224,6 +224,9 @@ ClassicalOrthogonalPolynomials.ShuffledFFT
 ClassicalOrthogonalPolynomials.legendre_grammatrix
 ```
 ```@docs
+ClassicalOrthogonalPolynomials.weightedgrammatrix
+```
+```@docs
 ClassicalOrthogonalPolynomials.WeightedOPLayout
 ```
 ```@docs

examples/pfem.jl

@@ -0,0 +1,67 @@
+using ClassicalOrthogonalPolynomials, Plots, Test
+
+
+#
+# We can solve ODEs like an inhomogenous Airy equation
+# $$
+# u'' - x*u = f
+# $$
+# with zero Dirichlet conditions using a basis built from ultraspherical polynomials
+# of the form
+# $$
+# (1-x^2)*C_n^{(3/2)}(x) = c_n (1-x^2) P_n^{(1,1)}(x)
+# $$
+# We construct this basis as follows:
+
+C = Ultraspherical(3/2)
+W = Weighted(C)
+plot(W[:,1:4])
+
+# We can differentiate using the `diff` command. Unlike arrays, for quasi-arrays `diff(W)` defaults to
+# `diff(W; dims=1)`. 
+
+plot(diff(W)[:,1:4])
+
+# We can get out a differentiation matrix via
+
+D² = C \ diff(diff(W))
+
+# We can construct the multiplication by $x$ with the connection between `W` and `C`
+# via:
+
+x = axes(W,1)
+X = C \ (x .* W)
+
+# Thus our operator becomes:
+
+L = D² - X
+
+# We can compute the coefficients using:
+
+c = L \ transform(C, exp)
+u = W*c
+plot(u)
+
+
+# Weak form for the Laplacian with W as a basis for test/trial  is given by
+# $$
+# ⟨dv/dx, du/dx ⟩ ≅ diff(v)'diff(u) = diff(W*d)'*diff(W*c) == d'*diff(W)'diff(W)*c
+# $$
+# where $v = W*d$ for some vector of coeficients d.  For $x$ we have
+# $$
+# ⟨v, x*u ⟩ ≅ v'(x .* u) == d'*W'(x .* W)*c
+# $$
+
+Δ  = -(diff(W)'diff(W)) # or weaklaplacian(W)
+x = axes(W,1)
+X = W' * (x .* W)
+L = Δ - X
+
+# We can solve an ODE via:
+
+c = L \ W'exp.(x)
+u = W*c
+plot!(u)
+
+
+# The two solvers match!

src/ClassicalOrthogonalPolynomials.jl

@@ -226,6 +226,9 @@ orthogonalityweight(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVe
 
 weighted(P::AbstractQuasiMatrix) = Weighted(P)
 
+"""
+gives the inner products of OPs with their weight, i.e., Weighted(P)'P.
+"""
 weightedgrammatrix(P) = weightedgrammatrix_layout(MemoryLayout(P), P)
 function weightedgrammatrix_layout(::MappedOPLayout, P)
     Q = parent(P)
@@ -292,6 +295,16 @@ plotgrid_layout(::AbstractOPLayout, P, n::Integer) = grid(P, min(40n, MAX_PLOT_P
 plotgrid_layout(::MappedOPLayout, P, n::Integer) = plotgrid_layout(MappedBasisLayout(), P, n)
 plotvalues_layout(::ExpansionLayout{MappedOPLayout}, f, x...) = plotvalues_layout(ExpansionLayout{MappedBasisLayout}(), f, x...)
 
+hasboundedendpoints(_) = false # assume blow up
+function plotgrid_layout(::WeightedOPLayout, P, n::Integer)
+    if hasboundedendpoints(weight(P))
+        plotgrid(unweighted(P), n)
+    else
+        grid(unweighted(P), min(40n, MAX_PLOT_POINTS))
+    end
+end
+
+
 function golubwelsch(X::AbstractMatrix)
     D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition
     D, V[1,:].^2

src/classical/chebyshev.jl

@@ -12,9 +12,13 @@ ChebyshevWeight() = ChebyshevWeight{1,Float64}()
 AbstractQuasiArray{T}(::ChebyshevWeight{kind}) where {T,kind} = ChebyshevWeight{kind,T}()
 AbstractQuasiVector{T}(::ChebyshevWeight{kind}) where {T,kind} = ChebyshevWeight{kind,T}()
 
+const ChebyshevTWeight = ChebyshevWeight{1}
+const ChebyshevUWeight = ChebyshevWeight{2}
+
+getproperty(w::ChebyshevTWeight{T}, ::Symbol) where T = -one(T)/2
+getproperty(w::ChebyshevUWeight{T}, ::Symbol) where T = one(T)/2
 
-getproperty(w::ChebyshevWeight{1,T}, ::Symbol) where T = -one(T)/2
-getproperty(w::ChebyshevWeight{2,T}, ::Symbol) where T = one(T)/2
+hasboundedendpoints(::ChebyshevUWeight) = true
 
 """
 Chebyshev{kind,T}()
@@ -28,8 +32,6 @@ Chebyshev{kind}() where kind = Chebyshev{kind,Float64}()
 AbstractQuasiArray{T}(::Chebyshev{kind}) where {T,kind} = Chebyshev{kind,T}()
 AbstractQuasiMatrix{T}(::Chebyshev{kind}) where {T,kind} = Chebyshev{kind,T}()
 
-const ChebyshevTWeight = ChebyshevWeight{1}
-const ChebyshevUWeight = ChebyshevWeight{2}
 const ChebyshevT = Chebyshev{1}
 const ChebyshevU = Chebyshev{2}
 

src/classical/jacobi.jl

@@ -46,6 +46,8 @@ summary(io::IO, w::JacobiWeight) = print(io, "(1-x)^$(w.a) * (1+x)^$(w.b) on -1.
 
 sum(P::JacobiWeight) = jacobimoment(P.a, P.b)
 
+hasboundedendpoints(w::AbstractJacobiWeight) = w.a ≥ 0 && w.b ≥ 0
+
 
 # support auto-basis determination
 

src/classical/ultraspherical.jl

@@ -28,6 +28,8 @@ end
 
 sum(w::UltrasphericalWeight{T}) where T = sqrt(convert(T,π))*exp(loggamma(one(T)/2 + w.λ)-loggamma(1+w.λ))
 
+hasboundedendpoints(w::UltrasphericalWeight) = 2w.λ ≥ 1
+
 
 struct Ultraspherical{T,Λ} <: AbstractJacobi{T}
     λ::Λ
@@ -80,6 +82,19 @@ plan_transform(P::Ultraspherical{T}, szs::NTuple{N,Int}, dims...) where {T,N} =
 Jacobi(C::Ultraspherical{T}) where T = Jacobi(C.λ-one(T)/2,C.λ-one(T)/2)
 
 
+
+######
+# Weighted Gram Matrix
+######
+
+# 2^(1-2λ)*π*gamma(n+2λ)/((n+λ)*gamma(λ)^2 * n!)
+function weightedgrammatrix(P::Ultraspherical{T}) where T
+    λ = P.λ
+    n = 0:∞
+    c = 2^(1-2λ) * convert(T,π)/gamma(λ)^2
+    Diagonal(c * exp.(loggamma.(n .+ 2λ) .- loggamma.(n .+ 1) ) ./ (n .+ λ))
+end
+
 ########
 # Jacobi Matrix
 ########
@@ -129,7 +144,11 @@ function diff(WS::Weighted{T,<:Ultraspherical}; dims=1) where T
     else
         n = (0:∞)
         A = _BandedMatrix((-one(T)/(2*(λ-1)) * ((n.+1) .* (n .+ (2λ-1))))', ℵ₀, 1,-1)
-        ApplyQuasiMatrix(*, Weighted(Ultraspherical{T}(λ-1)), A)
+        if λ == 3/2
+            ApplyQuasiMatrix(*, Legendre{T}(), A)
+        else
+            ApplyQuasiMatrix(*, Weighted(Ultraspherical{T}(λ-1)), A)
+        end
     end
 end
 
@@ -174,6 +193,8 @@ function \(U::Ultraspherical{<:Any,<:Integer}, C::ChebyshevU)
     U\Ultraspherical(C)
 end
 
+
+
 \(T::Chebyshev, C::Ultraspherical) = inv(C \ T)
 
 function \(C2::Ultraspherical{<:Any,<:Integer}, C1::Ultraspherical{<:Any,<:Integer})
@@ -209,29 +230,48 @@ function \(C2::Ultraspherical, C1::Ultraspherical)
     end
 end
 
-function \(w_A::WeightedUltraspherical, w_B::WeightedUltraspherical)
-    wA,A = w_A.args
-    wB,B = w_B.args
+function \(w_A::Weighted{<:Any,<:Ultraspherical}, w_B::Weighted{<:Any,<:Ultraspherical})
+    A = w_A.P
+    B = w_B.P
     T = promote_type(eltype(w_A),eltype(w_B))
 
-    if wA == wB
-        A \ B
-    elseif B.λ == A.λ+1 && wB.λ == wA.λ+1 # Lower
+    if A == B
+        SquareEye{T}(ℵ₀)
+    elseif B.λ == A.λ+1
         λ = convert(T,A.λ)
         _BandedMatrix(Vcat(((2λ:∞) .* ((2λ+1):∞) ./ (4λ .* (λ+1:∞)))',
                             Zeros{T}(1,∞),
                             (-(1:∞) .* (2:∞) ./ (4λ .* (λ+1:∞)))'), ℵ₀, 2,0)
+    elseif B.λ > A.λ+1
+        J = Weighted(Ultraspherical(B.λ-1))
+        (w_A\J) * (J\w_B)
     else
-        error("not implemented for $A and $wB")
+        error("not implemented for $w_A and $w_B")
     end
 end
 
-\(w_A::Weighted{<:Any,<:Ultraspherical}, w_B::Weighted{<:Any,<:Ultraspherical}) = convert(WeightedBasis, w_A) \ convert(WeightedBasis, w_B)
 
-\(A::Legendre, wB::WeightedUltraspherical) = Ultraspherical(A) \ wB
 
-function \(A::Ultraspherical, w_B::WeightedUltraspherical) 
-    (UltrasphericalWeight(zero(A.λ)) .* A) \ w_B
+function \(w_A::WeightedUltraspherical, w_B::WeightedUltraspherical)
+    wA,A = w_A.args
+    wB,B = w_B.args
+    T = promote_type(eltype(w_A),eltype(w_B))
+
+    if wA == wB
+        A \ B
+    elseif wA.λ == A.λ && wB.λ == B.λ # weighted
+        Weighted(A) \ Weighted(B)
+    elseif wB.λ ≥ wA.λ+1 # lower
+        J = UltrasphericalWeight(wB.λ-1) .* Ultraspherical(B.λ-1)
+        (w_A\J) * (J\w_B)
+    else
+        error("not implemented for $w_A and $w_B")
+    end
 end
 
+\(w_A::WeightedUltraspherical, w_B::Weighted{<:Any,<:Ultraspherical}) = w_A \ convert(WeightedBasis,w_B)
+\(w_A::Weighted{<:Any,<:Ultraspherical}, w_B::WeightedUltraspherical) = convert(WeightedBasis,w_A) \ w_B
+\(A::Ultraspherical, w_B::Weighted{<:Any,<:Ultraspherical}) = A \ convert(WeightedBasis,w_B)
+\(A::Ultraspherical, w_B::WeightedUltraspherical) = (UltrasphericalWeight(one(A.λ)/2) .* A) \ w_B
+\(A::Legendre, wB::WeightedUltraspherical) = Ultraspherical(A) \ wB
 

test/runtests.jl

@@ -78,4 +78,4 @@ end
     struct MyIncompleteJacobi <: ClassicalOrthogonalPolynomials.AbstractJacobi{Float64} end
     @test_throws ErrorException jacobimatrix(MyIncompleteJacobi())
     @test_throws ErrorException plan_transform(MyIncompleteJacobi(), 5)
-end
+end

test/test_odes.jl

@@ -1,9 +1,11 @@
 using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices,
-        SemiseparableMatrices, LazyArrays, ArrayLayouts, Test
+        LazyArrays, ArrayLayouts, Test
+
+using SemiseparableMatrices
 
 import QuasiArrays: MulQuasiMatrix
 import ClassicalOrthogonalPolynomials: oneto
-import ContinuumArrays: MappedBasisLayout, MappedWeightedBasisLayout
+import ContinuumArrays: MappedBasisLayout, MappedWeightedBasisLayout, weaklaplacian
 import LazyArrays: arguments, ApplyMatrix, colsupport, MemoryLayout
 import SemiseparableMatrices: VcatAlmostBandedLayout
 
@@ -246,4 +248,31 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
             @test u[0.1] ≈ 1.6878004187402804
         end
     end
+
+
+    @testset "ultraspherical ODE" begin
+        @testset "strong form" begin
+            C = Ultraspherical(3/2)
+            W = Weighted(C)
+            D² = C \ diff(diff(W))
+            x = axes(W,1)
+            X = C \ (x .* W)
+            L = D² - X
+            c = L \ transform(C, exp)
+            u = W*c
+            @test u[0.0] ≈ -0.536964648316337 # mathematica
+        end
+
+        @testset "weak form" begin
+            C = Ultraspherical(3/2)
+            W = Weighted(C)
+            Δ = weaklaplacian(W)
+            x = axes(W,1)
+            X = W' * (x .* W)
+            L = Δ - X
+            c = L \ (W'exp.(x))
+            u = W*c
+            @test u[0.0] ≈ -0.536964648316337 # mathematica
+        end
+    end
 end