Skip to content

Jacobi matrices for RectPolynomial #173

@MikaelSlevinsky

Description

@MikaelSlevinsky

I just added these in commit b2f95fe, but they cannot be instantiated like T² \ (x .* T²), like those on the Dunkl-Xu disk or the JacobiTriangle for reasons I don't understand. The code is so simple that the error must be somewhere else. I suspect it's because the Jacobi matrices for the basis are lazy tridiagonals instead of banded-block-banded matrices for the nontrivial Koornwinder classes, which means the error is not really in this code.

julia> using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials

julia> J = JacobiTriangle(0, 0, 0)
JacobiTriangle(0, 0, 0)

julia> jacobimatrix(Val{1}(), J); # works

julia> x,y = first.(axes(J,1)), last.(axes(J,1))
(first.(Inclusion(UnitSimplex(Val(2)))), last.(Inclusion(UnitSimplex(Val(2)))))

julia> jacobimatrix(Val{1}(), J)^C

julia> J \ (x .* J); # also works

julia>= RectPolynomial(ChebyshevT(), ChebyshevT())
ChebyshevT()  ChebyshevT()

julia> jacobimatrix(Val{1}(), T²); # works. This is what I just added

julia> x,y = first.(axes(T²,1)), last.(axes(T²,1))
(first.(Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev)))), last.(Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev)))))

julia>\ (x .* T²) # errs
ERROR: DimensionMismatch: Second axis of A, 1:1:ℵ₀, and first axis of B, Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev))) must match
Stacktrace:
 [1] _check_mul_axes(A::RectPolynomial{Float64, Tuple{…}}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:95
 [2] check_mul_axes(::RectPolynomial{Float64, Tuple{…}}, ::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:111
 [3] instantiate
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:123 [inlined]
 [4] materialize
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:127 [inlined]
 [5] mul
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:128 [inlined]
 [6] *(A::RectPolynomial{Float64, Tuple{…}}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/ZCn0F/src/matmul.jl:23
 [7] broadcasted(::QuasiArrays.LazyQuasiArrayStyle{…}, ::typeof(*), x::QuasiArrays.BroadcastQuasiVector{…}, P::RectPolynomial{…})
   @ HarmonicOrthogonalPolynomials ~/.julia/packages/HarmonicOrthogonalPolynomials/ViYwn/src/multivariateops.jl:45
 [8] broadcasted(::typeof(*), ::QuasiArrays.BroadcastQuasiVector{…}, ::RectPolynomial{…})
   @ Base.Broadcast ./broadcast.jl:1347
 [9] top-level scope
   @ REPL[119]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> 

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions