Ideal Fluid Flow off an interval (#43)

dlfivefifty · web-flow · commit d9a76016b798 · 2021-06-25T11:04:06.000+01:00
* Add IdealFluidFlow example

* Update test_stieltjes.jl

* tests fixed

* Add sum and hilbert special case

* Update Project.toml

* increase coverafge

* increase coverage

* increase coverage
diff --git a/Project.toml b/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"
diff --git a/examples/idealfluidflow.jl b/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))
diff --git a/examples/ultrasphericalspectralmethod.jl b/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)
diff --git a/src/ClassicalOrthogonalPolynomials.jl b/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)
diff --git a/src/classical/chebyshev.jl b/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}}
 
diff --git a/src/classical/jacobi.jl b/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
diff --git a/src/clenshaw.jl b/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)
 
diff --git a/src/normalized.jl b/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)
diff --git a/src/stieltjes.jl b/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)
diff --git a/test/test_chebyshev.jl b/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
diff --git a/test/test_normalized.jl b/test/test_normalized.jl
@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, ArrayLayouts, LazyArrays, Base64, LinearAlgebra, Test
-import ClassicalOrthogonalPolynomials: NormalizedBasisLayout, recurrencecoefficients, Normalized, Clenshaw, weighted
+import ClassicalOrthogonalPolynomials: NormalizedBasisLayout, recurrencecoefficients, Normalized, Clenshaw, weighted, grid, plotgrid
 import LazyArrays: CachedVector, PaddedLayout
 import ContinuumArrays: MappedWeightedBasisLayout
 
@@ -202,4 +202,11 @@ import ContinuumArrays: MappedWeightedBasisLayout
         # Q[x,:]' * ((X*Pn - Pn*X)* Q[y,:])
         @test (x-y) * Q[x,1:n]'*Q[y,1:n] ≈ Q[x,n:n+1]' * (X*Pn - Pn*X)[n:n+1,n:n+1] * Q[y,n:n+1]
     end
+
+    @testset "plotting" begin
+        P = Legendre()
+        Q = Normalized(P)
+        @test grid(Q[:,1:5]) == grid(Q[:,collect(1:5)]) == grid(P[:,1:5])
+        @test plotgrid(Q[:,1:5]) == plotgrid(Q[:,collect(1:5)]) == plotgrid(P[:,1:5])
+    end
 end
diff --git a/test/test_odes.jl b/test/test_odes.jl
@@ -94,7 +94,7 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
         S = Jacobi(true,true)
         w = JacobiWeight(true,true)
         D = Derivative(axes(S,1))
-        X = QuasiDiagonal(Inclusion(axes(S,1)))
+        X = QuasiDiagonal(axes(S,1))
 
         @test ((Legendre() \ S)*(S\(w.*S)))[1:10,1:10] ≈ (Legendre() \ (w .* S))[1:10,1:10]
         @test (Ultraspherical(3/2)\(D^2*(w.*S)))[1:10,1:10] ≈ diagm(0 => -(2:2:20))
@@ -169,6 +169,13 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
         L = Vcat(T[[-1,1],:], A)
         u = L \ [airyai(-ε^(-2/3)); airyai(ε^(2/3)); zeros(∞)]
         @test T[-0.1,:]'u ≈ airyai(-0.1*ε^(-2/3))
+
+        @testset "auto-vcat" begin
+            L = [T[[begin,end],:]; C \ (D^2 * T - x .* T)]
+            c = L \ [airyai(-1); airyai(1); zeros(∞)]
+            u = T*c
+            @test u[0.0] ≈ airyai(0.0)
+        end
     end
 
     @testset "combo operators" begin
diff --git a/test/test_stieltjes.jl b/test/test_stieltjes.jl
@@ -36,8 +36,8 @@ end
 
         f = wT * [[1,2,3]; zeros(∞)];
         J = T \ (x .* T)
-        # TODO: fix LazyBandedMatrices.copy(M::Mul{Broadcast...}) to call simplify
-        @test_broken π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f)[1]
+    
+        @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f)[1]
         @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f.args[1]*f.args[2])[1]
 
         x = Inclusion(0..1)
@@ -47,6 +47,19 @@ end
         f = wT2 * [[1,2,3]; zeros(∞)];
 
         @test (π/2*(((z-1/2)*I - J/2) \ f.args[2]))[1] ≈ (S*f.args[1]*f.args[2])[1]
+
+        @testset "Real point" begin
+            t = 2.0
+            T = ChebyshevT()
+            U = ChebyshevU()
+            x = axes(T,1)
+            @test inv.(t .- x') * Weighted(T) ≈ inv.((t+eps()im) .- x') * Weighted(T)
+            @test (inv.(t .- x') * Weighted(U))[1:10] ≈ (inv.((t+eps()im) .- x') * Weighted(U))[1:10]
+            
+            t = 0.5
+            @test_broken inv.(t .- x') * Weighted(T)
+            @test_broken inv.(t .- x') * Weighted(U)
+        end
     end
 
     @testset "Hilbert" begin
@@ -108,7 +121,7 @@ end
         L = log.(abs.(x .- x'))
         u =  wT * (2 *(T \ exp.(x)))
         @test u[0.1] ≈ exp(0.1)/sqrt(0.1-0.1^2)
-        @test_broken (L * u)[0.5] ≈ -7.471469928754152 # Mathematica
+        @test (L * u)[0.5] ≈ -7.471469928754152 # Mathematica
     end
 
     @testset "pow kernel" begin
@@ -118,6 +131,25 @@ end
         @test S isa ClassicalOrthogonalPolynomials.PowKernel
         @test_broken S*P
     end
+
+    @testset "Ideal Fluid Flow" begin
+        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(); 
+        @test (inv.(0.1+ε*im .- x') * u + inv.(0.1-ε*im .- x') * u)/2 ≈ imag(c*0.1)
+        @test 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))
+        @test v(0.1,0.2) ≈ 0.18496257285081724 # Emperical
+    end
 end
 
 #################################################
