Move over TwoBandJacobi (#82)

* Move over TwoBandJacobi

* Update Project.toml

* work on tests

* fix tests

* Update Project.toml

* Update Project.toml

* two band tests

* BandedMatrices v1

Project.toml

@@ -1,7 +1,7 @@
 name = "AlgebraicCurveOrthogonalPolynomials"
 uuid = "b2d6dd58-be48-4100-8375-7f7d22241ed8"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.0.1"
+version = "0.0.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -27,32 +27,34 @@ MultivariateOrthogonalPolynomials = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SemiclassicalOrthogonalPolynomials = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
+SingularIntegrals = "d7440221-8b5e-42fc-909c-0567823f424a"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-ArrayLayouts = "0.8, 1"
-BandedMatrices = "0.17"
+ArrayLayouts = "1.4.2"
+BandedMatrices = "0.17, 1"
 BlockArrays = "0.16"
-BlockBandedMatrices = "0.11, 0.12"
-ClassicalOrthogonalPolynomials = "0.7, 0.8, 0.9"
-ContinuumArrays = "0.12"
+BlockBandedMatrices = "0.12"
+ClassicalOrthogonalPolynomials = "0.11"
+ContinuumArrays = "0.16"
 DomainSets = "0.5, 0.6"
-FastGaussQuadrature = "0.4.3, 0.5"
-FastTransforms = "0.13, 0.14, 0.15"
+FastGaussQuadrature = "0.5"
+FastTransforms = "0.15"
 FillArrays = "0.13, 1"
 ForwardDiff = "0.10.12"
-HarmonicOrthogonalPolynomials = "0.4"
-InfiniteArrays = "0.12, 0.13"
-InfiniteLinearAlgebra = "0.6"
+HarmonicOrthogonalPolynomials = "0.5"
+InfiniteArrays = "0.13"
+InfiniteLinearAlgebra = "0.7.1"
 Infinities = "0.1"
-LazyArrays = "0.22.9, 1"
-LazyBandedMatrices = "0.8"
-MatrixFactorizations = "0.9, 1"
-MultivariateOrthogonalPolynomials = "0.4.1"
-QuasiArrays = "0.9, 0.10"
-SemiclassicalOrthogonalPolynomials = "0.3"
+LazyArrays = "1.8.1"
+LazyBandedMatrices = "0.9"
+MatrixFactorizations = "2"
+MultivariateOrthogonalPolynomials = "0.6"
+QuasiArrays = "0.11"
+SemiclassicalOrthogonalPolynomials = "0.5"
 SemiseparableMatrices = "0.3"
+SingularIntegrals = "0.2"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.7"

examples/equilibriummeasure.jl

@@ -0,0 +1,137 @@
+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))
+Base.big(d::Dual{T,V,N}) where {T,V,N} = Dual{T}(big(d.value), ForwardDiff.Partials(map(big,d.partials.values)))
+
+#### 
+#
+# 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
+
+μ = equilibrium(sqrt(2))
+
+T = Chebyshev()
+b = sqrt(2)
+μ = Weighted(T) * equilibriumcoefficients(T, b)
+x = axes(μ,1)
+
+plot(μ)
+
+xx = 0.7; 2b*(log.(abs.(x .- x'))*μ)[xx] - V(b*xx)
+
+
+
+b = 1.0 # initial guess
+for _ = 1:10
+    b -= derivative(equilibriumendpointvalue,b) \ equilibriumendpointvalue(b)
+end
+b
+
+plot(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)
+    x = axes(Q,1)
+    H = @. inv(x - t')
+    H̃ = Q \ H*W
+    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \( Q \ derivative.(V, b*x))]
+end
+function equilibriumendpointvalues(ab::SVector{2})
+    a,b = ab
+    # orthogonal polynomials w.r.t.
+    # abs(x) / (sqrt(1-x^2) * sqrt(x^2 - ρ^2))
+    P = TwoBandJacobi(a/b, -one(a)/2, -one(a)/2, one(a)/2)
+    Vector(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
+
+ab = SVector(2.,3.)
+ab -= jacobian(equilibriumendpointvalues,ab) \ equilibriumendpointvalues(ab)
+a,b = ab
+xx = range(-4,4;length=1000)
+μ = equilibrium(ab)
+μx = x -> a < abs(x) < b ? μ[x/b] : 0.0
+plot!(xx, μx.(xx))
+
+plot(equilibrium(ab))
+

examples/idealfluidflow_2_plates.jl

@@ -0,0 +1,17 @@
+using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, LinearAlgebra
+import ClassicalOrthogonalPolynomials: associated
+
+ρ = 0.5
+T = TwoBandJacobi(ρ,-1/2,-1/2,1/2)
+U = associated(T)
+x = axes(T,1)
+H = inv.(x .- x')
+
+L = U \ (H * Weighted(T))
+c = L[:,2:∞] \ (U \ x)
+c[1] = -c[2]
+
+u = Weighted(T) * c
+u[0.5000001]
+
+[0.1]

src/AlgebraicCurveOrthogonalPolynomials.jl

@@ -6,21 +6,24 @@ using FastGaussQuadrature, FastTransforms, SpecialFunctions, LinearAlgebra, Bloc
 import ForwardDiff: jacobian
 import ForwardDiff: jacobian, Dual, gradient, value, partials
 import FillArrays: SquareEye
-import LinearAlgebra: eigvals, eigen, isapprox, SymTridiagonal, norm, factorize
+import LinearAlgebra: eigvals, eigen, isapprox, norm, factorize
 import FastGaussQuadrature: jacobimoment
 import QuasiArrays: DefaultQuasiArrayStyle, cardinality, LazyQuasiArrayStyle
-import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff
-import ContinuumArrays: Weight, grid, ℵ₁, ℵ₀, @simplify, ProjectionFactorization, plan_grid_transform, unweighted, weight
-import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace
+import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff, copy, sum, size, axes
+import ContinuumArrays: Weight, grid, ℵ₁, ℵ₀, @simplify, ProjectionFactorization, plan_grid_transform, unweighted, weight, transform_ldiv
+import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace, WeightedBasis, HalfWeighted, plotgrid, golubwelsch, recurrencecoefficients
 import MultivariateOrthogonalPolynomials: BlockOneTo, ModalInterlace, BlockRange1, Plan, ModalTrav
 import BlockArrays: block, blockindex, _BlockedUnitRange, BlockSlice, blockcolsupport
-import BlockBandedMatrices: BlockTridiagonal, AbstractBlockBandedMatrix, blockbandwidths, subblockbandwidths, _BandedBlockBandedMatrix
-import SemiclassicalOrthogonalPolynomials: divmul, HalfWeighted
+import BlockBandedMatrices: BlockTridiagonal, AbstractBlockBandedMatrix, blockbandwidths, subblockbandwidths, _BandedBlockBandedMatrix, band
+import SemiclassicalOrthogonalPolynomials: HalfWeighted
+import LazyArrays: LazyVector
+import SingularIntegrals: Associated, associated, stieltjes
+import LazyBandedMatrices: Tridiagonal, SymTridiagonal, BandedMatrix
 
 export quarticjacobi, blocksymtricirculant, unroll, randspeccurve, speccurve, specgrid, speccurvemat, symroll, symunroll, spec2alg,
         wedgep, wedgeq, wedger, wedgetransform, plan_wedgetransform, plan_squaretransform, gausswedge, JacobiWedge, LegendreSquare, gausssquare,
         HermLaurent, ImHermLaurent, jointeigen, jointeigvals, BlockTridiagonal, LegendreCircle, UltrasphericalCircle, Block, SVector, CircleCoordinate,
-        UltrasphericalArc, LegendreCubic, ZernikeAnnulus, ComplexZernikeAnnulus, hermlaurent
+        UltrasphericalArc, LegendreCubic, hermlaurent, TwoBandWeight, TwoBandJacobi
 
 
 
@@ -103,10 +106,10 @@ include("wedge.jl")
 include("square.jl")
 include("cubic.jl")
 include("quartic.jl")
-include("annulus.jl")
 
 include("algcurvapprox.jl")
 
 include("hermlaurent.jl")
+include("twoband.jl")
 
 end # module