Ideal Fluid Flow off an interval (#43)

* Add IdealFluidFlow example

* Update test_stieltjes.jl

* tests fixed

* Add sum and hilbert special case

* Update Project.toml

* increase coverafge

* increase coverage

* increase coverage

Author: dlfivefifty

Project.toml

@@ -22,6 +22,7 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
@@ -36,11 +37,11 @@ FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.12"
 FillArrays = "0.11"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.10.3"
-InfiniteLinearAlgebra = "0.5.3"
+InfiniteArrays = "0.11.1"
+InfiniteLinearAlgebra = "0.5.8"
 IntervalSets = "0.5"
 LazyArrays = "0.21"
-LazyBandedMatrices = "0.5.5"
+LazyBandedMatrices = "0.6"
 QuasiArrays = "0.6"
 SpecialFunctions = "1.0"
 julia = "1.6"

examples/idealfluidflow.jl

@@ -0,0 +1,32 @@
+using ClassicalOrthogonalPolynomials, Plots
+
+##
+#  Ideal fluid flow consists of level sets of the imagainary part of a function
+# that is asymptotic to c*z and whose imaginary part vanishes on Γ
+#
+#
+# On the unit interval, -2*H*u gives the imaginary part of cauchy(u,x)
+#  So if we want to find u defined on Γ so that hilbert(u,x) = imag(c*x)
+#  then c*z + 2cauchy(u,z) vanishes on Γ
+##
+
+T = ChebyshevT()
+U = ChebyshevU()
+x = axes(U,1)
+H = inv.(x .- x')
+
+c = exp(0.5im)
+
+
+u = Weighted(U) * ((H * Weighted(U)) \ imag(c * x))
+
+ε  = eps(); (inv.(0.1+ε*im .- x') * u + inv.(0.1-ε*im .- x') * u)/2 ≈ imag(c*0.1)
+ε  = eps(); real(inv.(0.1+ε*im .- x') * u ) ≈ imag(c*0.1)
+
+v = (s,t) -> (z = (s + im*t); imag(c*z) - real(inv.(z .- x') * u))
+
+
+xx = range(-3,3; length=100)
+yy = range(-1,1; length=100)
+plot([-1,1],[0,0]; color=:black)
+contour!(xx,yy,v.(xx',yy))

examples/ultrasphericalspectralmethod.jl

@@ -1,11 +1,21 @@
-using ClassicalOrthogonalPolynomials, Plots
+using ClassicalOrthogonalPolynomials, Plots, Test
+
+###
+# We can solve ODEs like the Airy equation
+#
+# u(-1) = airyai(-1)
+# u(1) = airyai(1)
+# u'' = x * u
+#
+# using the ultraspherical spectral method. 
 
 T = Chebyshev()
 C = Ultraspherical(2)
-D = Derivative(axes(T,1))
-A = (C\(D^2*T))+100(C\T)
+x = axes(T,1)
+D = Derivative(x)
 
-c = Vcat(T[[-1,1],:], A) \ [1;2;zeros(∞)]
+c = [T[[begin,end],:]; C \ (D^2 * T - x .* T)] \ [airyai(-1); airyai(1); zeros(∞)]
 u = T*c
 
+@test u[0.0] ≈ airyai(0.0)
 plot(u)

src/ClassicalOrthogonalPolynomials.jl

@@ -25,7 +25,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
-                    _getindex, layout_getindex, _factorize
+                    _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
@@ -59,7 +59,8 @@ include("interlace.jl")
 cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
 cardinality(::EuclideanDomain) = ℵ₁
 
-transform_ldiv(A, f, ::Tuple{<:Any,InfiniteCardinal{0}})  = adaptivetransform_ldiv(A, f)
+transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  = 
+    adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)
 
 function chop!(c::AbstractVector, tol::Real)
     @assert tol >= 0
@@ -279,7 +280,7 @@ _tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)
 jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) =
     _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))
 
-grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,AbstractUnitRange}}) =
+grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}}) =
     eigvals(symtridiagonalize(jacobimatrix(P)))
 
 function golubwelsch(X)

src/classical/chebyshev.jl

@@ -22,6 +22,8 @@ on -1..1.
 struct Chebyshev{kind,T} <: AbstractJacobi{T} end
 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 WeightedChebyshev{kind,T} = WeightedBasis{T,<:ChebyshevWeight{kind},<:Chebyshev{kind}}
 

src/classical/jacobi.jl

@@ -204,12 +204,12 @@ summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
 # transforms
 ###
 
-function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,AbstractUnitRange}}) where T
+function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,Any}}) where T
     kr,jr = parentindices(Pn)
     ChebyshevGrid{1,T}(maximum(jr))
 end
 
-function plotgrid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,AbstractUnitRange}}) where T
+function plotgrid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{Inclusion,Any}}) where T
     kr,jr = parentindices(Pn)
     ChebyshevGrid{2,T}(40maximum(jr))
 end

src/clenshaw.jl

@@ -80,6 +80,7 @@ Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, ::Colon) = Base.unsafe_
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::Number) = Base.unsafe_getindex(P,x,oneto(n))[end]
 
 getindex(P::OrthogonalPolynomial, x::Number, jr::AbstractInfUnitRange{Int}) = view(P, x, jr)
+getindex(P::OrthogonalPolynomial, x::AbstractVector, jr::AbstractInfUnitRange{Int}) = view(P, x, jr)
 Base.unsafe_getindex(P::OrthogonalPolynomial{T}, x::Number, jr::AbstractInfUnitRange{Int}) where T = 
     BroadcastVector{T}(Base.unsafe_getindex, Ref(P), x, jr)
 

src/normalized.jl

@@ -83,7 +83,8 @@ QuasiArrays.ApplyQuasiArray(Q::Normalized) = ApplyQuasiArray(*, arguments(ApplyL
 
 ArrayLayouts.mul(Q::Normalized, C::AbstractArray) = ApplyQuasiArray(*, Q, C)
 
-grid(Q::SubQuasiArray{<:Any,2,<:Normalized}) = grid(view(parent(Q).P, parentindices(Q)...))
+grid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = grid(view(parent(Q).P, parentindices(Q)...))
+plotgrid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = plotgrid(view(parent(Q).P, parentindices(Q)...))
 
 # transform_ldiv(Q::Normalized, C::AbstractQuasiArray) = Q.scaling .\ (Q.P \ C)
 function transform_ldiv(Q::Normalized, C::AbstractQuasiArray)
@@ -215,3 +216,10 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::Weighted) =
 summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")
 
 __sum(::NormalizedBasisLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
+function __sum(::WeightedOPLayout, A, dims)
+    @assert dims == 1
+    Hcat(sum(weight(A)), Zeros{eltype(A)}(1,∞))
+end
+
+_sum(p::SubQuasiArray{T,1,<:Weighted,<:Tuple{Inclusion,Int}}, ::Colon) where T = 
+    parentindices(p)[2] == 1 ? convert(T, sum(weight(parent(p)))) : zero(T)

src/stieltjes.jl

@@ -143,9 +143,20 @@ end
     w = orthogonalityweight(P)
     X = jacobimatrix(P)
     z, x = parent(S).args[1].args
+    z in axes(P,1) && transpose((inv.(x .- x') * wP)[z,:])
     transpose((X'-z*I) \ [-sum(w)*_p0(P); zeros(∞)])
 end
 
+sqrtx2(z::Number) = sqrt(z-1)*sqrt(z+1)
+sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)
+
+@simplify function *(S::StieltjesPoint, wP::Weighted{<:Any,<:ChebyshevU})
+    z, x = parent(S).args[1].args
+    z in axes(wP,1) && transpose((inv.(x .- x') * wP)[z,:])
+    ξ = inv(z + sqrtx2(z))
+    transpose(π * ξ.^oneto(∞))
+end
+
 
 @simplify function *(S::StieltjesPoint, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(wT)

test/test_chebyshev.jl

@@ -32,6 +32,10 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @test axes(T[1:1,:]) === (oneto(1), oneto(∞))
             @test T[1:1,:][:,1:5] == ones(1,5)
             @test T[0.1,:][1:10] ≈ T[0.1,1:10] ≈ (T')[1:10,0.1]
+
+            @testset "inf-range-indexing" begin
+                @test T[[begin,end],2:∞][:,2:5] == T[[-1,1],3:6]
+            end
         end
 
         @testset "U" begin
@@ -136,6 +140,9 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @test WT \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ wT \ (exp.(x) ./ sqrt.(1 .- x.^2))
             @test WT[:,1:20] \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ (WT \ (exp.(x) ./ sqrt.(1 .- x.^2)))[1:20]
             @test WT \ (x .* WT) == T \ (x .* T)
+            @test sum(WT; dims=1)[:,1:10] ≈ [π zeros(1,9)]
+            @test sum(WT[:,1]) ≈ π
+            @test iszero(sum(WT[:,2]))
         end
 
         @testset "mapped" begin
@@ -359,6 +366,11 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
     @testset "plot" begin
         @test ContinuumArrays.plotgrid(ChebyshevT()[:,1:5]) == ChebyshevGrid{2}(200)
     end
+
+    @testset "conversion" begin
+        T = ChebyshevT()
+        @test ChebyshevT{ComplexF64}() ≡ convert(AbstractQuasiArray{ComplexF64}, T) ≡ convert(AbstractQuasiMatrix{ComplexF64}, T)
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end