Work on two-interval equilibrium measure (#42)

* Work on two-interval equilibrium measure

* two-band equilibrium measure converges

* update versions

* fix tests

* Increase coverage

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '^1.6.0-0'
+          - '1.6'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,14 +1,15 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 ContinuumArrays = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
@@ -17,17 +18,17 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.6, 0.7"
+ArrayLayouts = "0.7"
 BandedMatrices = "0.16"
-ClassicalOrthogonalPolynomials = "0.3.5, 0.4"
-ContinuumArrays = "0.7, 0.8"
+ClassicalOrthogonalPolynomials = "0.4.2"
+ContinuumArrays = "0.8"
 FillArrays = "0.11.5"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.10.2, 0.11"
+InfiniteArrays = "0.11"
 LazyArrays = "0.21.3"
-QuasiArrays = "0.5, 0.6"
-SpecialFunctions = "0.10, 1.0"
-julia = "1.5"
+QuasiArrays = "0.6"
+SpecialFunctions = "1.0"
+julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

examples/equilibriummeasure.jl

@@ -0,0 +1,117 @@
+using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, ForwardDiff, Plots, StaticArrays
+import ForwardDiff: derivative, jacobian, Dual
+import SemiclassicalOrthogonalPolynomials: Weighted
+import ClassicalOrthogonalPolynomials: associated, affine
+
+Base.floatmin(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,V,N}(floatmin(V))
+
+
+####
+#
+# We first do a single interval.
+# Equilibrium measures for a symmetric potential
+#  with one interval of support 
+# a measure w(x) supported on
+# [-b,b]
+# such that
+#       1. H*w == V'
+#       2. sum(w) == 1
+#       3.  w is bounded
+# 
+#  rescaling x == b*t these becomes find 
+# a measure w̃(t) supported on
+# [-1,1]
+# such that
+#       1. H*w̃ == V'(b*t)
+#       2. sum(w̃) == 1/b
+#       3.  w is bounded
+#
+# Note (1) and (2) can always be satisfied
+# thus the constraint comes from (3).
+# The following gives the evaluation of the
+# unweighted-component of the measure evaluated
+# at (a,b)
+#####
+
+
+V = x -> x^2
+function equilibriumcoefficients(T, b)
+    U = associated(T)
+    W = Weighted(T)
+    t = axes(W,1)
+    H = @. inv(t - t')
+    H̃ = U \ (H*W)
+    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \ ( U \ derivative.(V, b*t))]
+end
+function equilibriumendpointvalue(b::Number)
+    T = ChebyshevT{typeof(b)}()
+    dot(T[end,:], equilibriumcoefficients(T,b))
+end
+
+function equilibrium(b::Number)
+    T = ChebyshevT{typeof(b)}()
+    U = ChebyshevU{typeof(b)}()
+    # convert to Weighted(U) to make value at ±b accurate
+    Weighted(U)[affine(-b..b,axes(T,1)),:] * ((Weighted(T) \ Weighted(U))[3:end,:] \ equilibriumcoefficients(T,b)[3:end])
+end
+
+b = 1.0 # initial guess
+for _ = 1:10
+    b -= derivative(equilibriumendpointvalue,b) \ equilibriumendpointvalue(b)
+end
+
+plot(equilibrium(b))
+ClassicalOrthogonalPolynomials.plotgrid(equilibrium(b))
+
+
+#####
+# Equilibrium measures for a symmetric potential 
+# with two intervals of support consists of finding 
+# a measure w(x) supported on
+# [-b,-a] ∪ [a,b]
+# such that
+#       1. H*w == V'
+#       2. sum(w) == 1
+#       3.  w is bounded
+# 
+#  rescaling x == b*t these becomes find 
+# a measure w̃(t) supported on
+# [-1,-a/b] ∪ [a/b,1]
+# such that
+#       1. H*w̃ == V'(b*t)
+#       2. sum(w̃) == 1/b
+#       3.  w is bounded
+#
+# Note (1) and (2) can always be satisfied
+# thus the two constraints come from (3).
+# The following gives the evaluation of the
+# unweighted-component of the measure evaluated
+# at (a,b)
+#####
+V = x -> x^4 - 10x^2
+function equilibriumcoefficients(P,a,b)
+    W = Weighted(P)
+    Q = associated(P)
+    t = axes(W,1)
+    H = @. inv(t - t')
+    H̃ = Q \ (H*W)
+    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \( Q \ derivative.(V, b*t))]
+end
+function equilibriumendpointvalues(ab::SVector{2})
+    a,b = ab
+    P = TwoBandJacobi(a/b, -1/2, -1/2, 1/2)
+    P[[a/b,1],:] * equilibriumcoefficients(P,a,b)
+end
+
+function equilibrium(ab)
+    a,b = ab
+    P = TwoBandJacobi(a/b, -1/2, -1/2, 1/2)
+    Weighted(P) * equilibriumcoefficients(P,a,b)
+end
+
+(a,b) = ab = SVector(2.,3.)
+ab -= jacobian(equilibriumendpointvalues,ab) \ equilibriumendpointvalues(ab)
+
+
+plot(equilibrium(ab))
+

src/SemiclassicalOrthogonalPolynomials.jl

@@ -1,4 +1,5 @@
 module SemiclassicalOrthogonalPolynomials
+using ClassicalOrthogonalPolynomials: WeightedOPLayout
 using ClassicalOrthogonalPolynomials, FillArrays, LazyArrays, ArrayLayouts, QuasiArrays, InfiniteArrays, ContinuumArrays, LinearAlgebra, BandedMatrices,
         SpecialFunctions, HypergeometricFunctions
 
@@ -12,7 +13,7 @@ import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficie
                                         OrthogonalPolynomialRatio, Weighted, Expansion, UnionDomain, oneto, Hilbert, WeightedOrthogonalPolynomial, HalfWeighted,
                                         Associated
 import InfiniteArrays: OneToInf, InfUnitRange
-import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid
+import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid
 import FillArrays: SquareEye
 import HypergeometricFunctions: _₂F₁general2
 
@@ -297,13 +298,16 @@ copy(L::Ldiv{SemiclassicalJacobiLayout,<:NormalizedBasisLayout}) = copy(Ldiv{Sem
 copy(L::Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},BasisLayout}(L.A, L.B))
 copy(L::Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
 
-
 copy(L::Ldiv{MappedBasisLayout,SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
 copy(L::Ldiv{SemiclassicalJacobiLayout,MappedBasisLayout}) = semijacobi_ldiv(L.A, L.B)
 
+copy(L::Ldiv{WeightedOPLayout,SemiclassicalJacobiLayout}) = copy(Ldiv{WeightedOPLayout,BasisLayout}(L.A, L.B))
+copy(L::Ldiv{SemiclassicalJacobiLayout,WeightedOPLayout}) = copy(Ldiv{BasisLayout,WeightedOPLayout}(L.A, L.B))
+
 
 copy(L::Ldiv{SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
 copy(L::Ldiv{<:Any,SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
+copy(L::Ldiv{<:AbstractBasisLayout,SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
 function copy(L::Ldiv{SemiclassicalJacobiLayout,SemiclassicalJacobiLayout})
     Q,P = L.A,L.B
     @assert Q.t == P.t

src/twoband.jl

@@ -4,7 +4,7 @@ twoband(ρ) = UnionDomain(-1..(-ρ), ρ..1)
 """
     TwoBandWeight(ρ, a, b, c)
 
-is a quasi-vector representing `|x|^(2c) * (x^2-ρ)^b * (1-x^2)^a`  for `ρ ≤ |x| ≤ 1`
+is a quasi-vector representing `|x|^(2c) * (x^2-ρ^2)^b * (1-x^2)^a`  for `ρ ≤ |x| ≤ 1`
 """
 struct TwoBandWeight{T} <: Weight{T}
     ρ::T
@@ -28,7 +28,7 @@ end
 
 function sum(w::TwoBandWeight{T}) where T
     if 2w.a == 2w.b == -1 && 2w.c == 1
-        convert(T,π)/2
+        convert(T,π)
     else
         error("not implemented")
 #         1/2 \[Pi] (-((\[Rho]^(1 + 2 b + 2 c)
@@ -112,7 +112,10 @@ function grid(L::SubQuasiArray{T,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{I
     [-g; g]
 end
 
-    
+function plotgrid(L::SubQuasiArray{T,2,<:TwoBandJacobi,<:Tuple{Inclusion,AbstractUnitRange}}) where T
+    g = plotgrid(legendre(parent(L).ρ .. 1)[:,parentindices(L)[2]])
+    [-g; g]
+end
 
 LinearAlgebra.factorize(L::SubQuasiArray{<:Any,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{Inclusion,OneTo}}) =
     ContinuumArrays._factorize(BasisLayout(), L)

test/runtests.jl

@@ -1,5 +1,5 @@
 using SemiclassicalOrthogonalPolynomials
-using ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, QuasiArrays, Test, LazyArrays, LinearAlgebra
+using ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, QuasiArrays, Test, LazyArrays, FillArrays, LinearAlgebra
 import BandedMatrices: _BandedMatrix
 import SemiclassicalOrthogonalPolynomials: op_lowering, RaisedOP
 import ClassicalOrthogonalPolynomials: recurrencecoefficients, orthogonalityweight, symtridiagonalize, Expansion
@@ -310,6 +310,12 @@ end
         @test (Q * (Q \ exp.(x)))[0.1] ≈ exp(0.1)
         @test (R * (R \ exp.(x)))[0.1] ≈ exp(0.1)
     end
+
+    @testset "conversion" begin
+        D = legendre(0..1) \ P
+        @test (P \ legendre(0..1))[1:10,1:10] == inv(Matrix(D[1:10,1:10]))
+        @test_broken (P \ Weighted(P)) == (Weighted(P) \ P) == Eye(∞)
+    end
 end
 
 @testset "Hierarchy" begin

test/test_twoband.jl

@@ -1,5 +1,5 @@
 using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test
-import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
+import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated, plotgrid
 
 @testset "Two Band" begin
     @testset "TwoBandWeight" begin
@@ -33,7 +33,7 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
         x = axes(w,1)
         H = inv.(x .- x')
         @test iszero(H*w)
-        @test sum(w) == π/2
+        @test sum(w) ≈ π
 
         T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
         Q = associated(T)
@@ -45,4 +45,10 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
 
         @test_broken H*TwoBandWeight(ρ, 1/2, 1/2, -1/2)
     end
+
+    @testset "plotgrid" begin
+        ρ = 0.5
+        P = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
+        @test all(x -> ρ ≤ abs(x) ≤ 1, plotgrid(P[:,1:5]))
+    end
 end