Support for orthogonal polynomial-based spaces in ApproxFun
ApproxFun.jl is a package for approximating functions. This package implements orthogonal polynomial based spaces to ApproxFun, including Chebyshev, Jacobi, Laguerre, and Hermite.
Support for orthogonal polynomial-based spaces in ApproxFun
License: MIT License
Support for orthogonal polynomial-based spaces in ApproxFun
ApproxFun.jl is a package for approximating functions. This package implements orthogonal polynomial based spaces to ApproxFun, including Chebyshev, Jacobi, Laguerre, and Hermite.
This test failure showed up in a log:
Check resizing: Test Failed at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:226
Expression: norm(((Laplacian() * u + 100u) - (A * u)[2]).coefficients) < 1.0e-9
Evaluated: 1.3708085290831766e-9 < 1.0e-9
Stacktrace:
[1] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:226
[2] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
[3] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:176
[4] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
[5] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:5
I don't think it is a real error so it might make sense to loosen the tolerance a bit.
julia> size(Multiplication(Fun(), Chebyshev()))
(ℵ₀, ℵ₀)
julia> size(Multiplication(Fun(), Legendre()))
(ℵ₀, ℵ₀)
julia> size(Multiplication(Fun(), Ultraspherical(1)))
(∞, ∞)
julia> size(Multiplication(Fun(), Ultraspherical(2)))
(ℵ₀, ℵ₀)
This seems to happen because the Toeplitz-Hankel operators are defined on Sequence spaces, and
julia> ApproxFunBase.dimension(SequenceSpace())
∞
However, this seems to be an implementation detail. Should size
use the dimension of the domain/range spaces instead?
julia> ApproxFunBase.dimension(Ultraspherical(1))
ℵ₀
It's possible to represent functions with infinite coefficients using InfiniteArrays.jl, however, evaluation is broken:
julia> f = Fun(Chebyshev(), -2.0 .^(-(1:∞)))
Fun(Chebyshev(),[-0.5, -0.25, -0.125, -0.0625, -0.03125, -0.015625, -0.0078125, -0.00390625, -0.00195313, -0.000976563 … ])
julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
[1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
[2] chebyshev_clenshaw(::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
[3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
[4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
[5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
[6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},Float64,BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
[7] top-level scope at none:0
It's actually a bit simpler to think about in the interval arithmetic setting where we only need a bound on evaluation (otherwise we would need to talk about tolerances). That is, in the format that #2 would return:
julia> f = Fun(Chebyshev(), BroadcastArray(IntervalArithmetic.Interval, -2.0 .^(-(1:∞)), 2.0 .^(-(1:∞))))
Fun(Chebyshev(),IntervalArithmetic.Interval{Float64}[[-0.5, 0.5], [-0.25, 0.25], [-0.125, 0.125], [-0.0625, 0.0625], [-0.03125, 0.03125], [-0.015625, 0.015625], [-0.0078125, 0.0078125], [-0.00390625, 0.00390625], [-0.00195313, 0.00195313], [-0.000976563, 0.000976563] … ])
julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
[1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
[2] chebyshev_clenshaw(::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
[3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
[4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
[5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
[6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},IntervalArithmetic.Interval{Float64},BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
[7] top-level scope at none:0
For Chebyshev it should be easy to work out bounds and override chebyshev_clenshaw(::BroadcastArray{...}, x)
to return these bounds. We could also then override chebyshev_clenshaw(a::Vcat, x)
to do Clenshaw recursively, using the result of chebyshev_clenshaw(a.arrays[end], x)
to determine the initial condition for the rest of Clenshaw. This would give rigorous evaluation of functions.
We should have as goal a 10-minute test suite per job on CI, per JuliaApproximation package. The tests in the ApproxFun ecosystem evolved from
The following would ideally work:
julia> using ApproxFun
julia> import IntervalArithmetic: @interval
julia> Base.sinpi(x::IntervalArithmetic.Interval) = sin(π*x)
julia> Fun(exp, @interval(0)..@interval(1))
ERROR: MethodError: no method matching plan_chebyshevtransform(::Array{IntervalArithmetic.Interval{Float64},1})
Closest candidates are:
plan_chebyshevtransform(::Array{D<:DualNumbers.Dual,1}; kind) where D<:DualNumbers.Dual at /Users/solver/Projects/ApproxFun.jl/src/Extras/dualnumbers.jl:35
plan_chebyshevtransform(::AbstractArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},1}; kind) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/solver/.julia/packages/FastTransforms/vEjxF/src/chebyshevtransform.jl:29
plan_chebyshevtransform(::AbstractArray{T<:Union{Complex{BigFloat}, BigFloat},1}; kind) where T<:Union{Complex{BigFloat}, BigFloat} at /Users/solver/Projects/ApproxFun.jl/src/Extras/fftGeneric.jl:49
Stacktrace:
[1] plan_transform(::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:65
[2] transform(::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}) at /Users/solver/Projects/ApproxFunBase.jl/src/Space.jl:427
[3] default_Fun(::Type{IntervalArithmetic.Interval{Float64}}, ::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}, ::Type{Val{false}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:44
[4] default_Fun(::ApproxFunBase.DFunction, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Int64, ::Type{Val{false}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:57
[5] default_Fun(::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Int64) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:72
[6] default_Fun(::ApproxFunBase.DFunction, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:121
[7] Fun(::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:177
[8] Fun(::Function, ::Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:173
[9] top-level scope at none:0
Ideally, the transform would return a data structure representing an infinite vector with exponential decay. It's actually already possible to construct such a thing using InfiniteArrays.jl:
julia> BroadcastArray(IntervalArithmetic.Interval, -2.0 .^(-(1:∞)), 2.0 .^(-(1:∞)))
BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}} with indices OneToInf():
[-0.5, 0.5]
[-0.25, 0.25]
[-0.125, 0.125]
[-0.0625, 0.0625]
[-0.03125, 0.03125]
[-0.015625, 0.015625]
[-0.0078125, 0.0078125]
[-0.00390625, 0.00390625]
[-0.00195313, 0.00195313]
[-0.000976563, 0.000976563]
[-0.000488282, 0.000488282]
[-0.000244141, 0.000244141]
[-0.000122071, 0.000122071]
[-6.10352e-05, 6.10352e-05]
[-3.05176e-05, 3.05176e-05]
[-1.52588e-05, 1.52588e-05]
[-7.6294e-06, 7.6294e-06]
[-3.8147e-06, 3.8147e-06]
⋮
This would be combined with a finite vector of coefficients using Vcat
. Some mathematical thought is needed on how to calculate this automatically, which would require bounding in a Bernstein ellipse.
@dpsanders FYI.
using ApproxFunOrthogonalPolynomials
# For reference
chebyshevpoints(Float64, 3, Val(1)) # [0.866, 0.0, -0.866]
chebyshevpoints(Float64, 3, Val(2)) # [1.0, 0.0, -1.0]
# These should be points corresponding to the Chebyshev polynomials of the second kind, I thought?
points(Ultraspherical(1, -1..1), 3) # [0.866, 0.0, -0.866]
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
Hi there, I was going through some examples and trying to learn how to work with ApproxFun when I ran into this.
I am not sure if this is a superficial issue of invalid inputs not being caught, or if there is some deeper problems, but here are my observations:
Note: The issue is the Conversion()
method, but I use an Operator
with →
here because it results in an immediately evaluated representation. The following issues occurred with all the Operator
s I tried.
I presume the following should not be allowed to evaluate
julia> Operator(Fun()) : Ultraspherical(2) → Ultraspherical(1)
TimesOperator : Ultraspherical(2) → Ultraspherical(1)
NaN NaN 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱ ⋱
and a stack overflow error is not ideal for the following
julia> Operator(Fun()) : Ultraspherical(3) → Ultraspherical(1.9)
ERROR: StackOverflowError:
Stacktrace:
[1] Conversion(::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunOrthogonalPolynomials\jojxS\src\Spaces\Ultraspherical\UltrasphericalOperators.jl:141
[2] Conversion(::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunOrthogonalPolynomials\jojxS\src\Spaces\Ultraspherical\UltrasphericalOperators.jl:148 (repeats 14488 times)
[3] promoterangespace(::ApproxFunBase.ConcreteMultiplication{Chebyshev{ChebyshevInterval{Float64},Float64},Ultraspherical{Int64,ChebyshevInterval{Float64},Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Int64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:83
[4] promoterangespace at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:79 [inlined]
[5] →(::ApproxFunBase.ConcreteMultiplication{Chebyshev{ChebyshevInterval{Float64},Float64},Ultraspherical{Int64,ChebyshevInterval{Float64},Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:74
[6] top-level scope at none:0
If we want to catch these two cases, I believe changing the elseif
s as in the following should be sufficient?
function Conversion(A::Ultraspherical,B::Ultraspherical)
a=order(A); b=order(B)
if b==a
ConversionWrapper(Operator(I,A))
# elseif a<b≤a+1 || b<a≤b+1
elseif -1≤b-a≤1 && b≠1
ApproxFun.ConcreteConversion(A,B)
# elseif b ≠ 1
elseif b-a>1
d=domain(A)
ConversionWrapper(TimesOperator(Conversion(Ultraspherical(b-1,d),B),
Conversion(A,Ultraspherical(b-1,d))))
else
error("Cannot convert from $A to $B")
end
end
I also noted some inconsistent behaviour of getindex()
when passing Ultrasphericals of Float and/or Rational order (inconsistency of when it throws an error, not inconsistencies in evaluation), but I suppose that is a separate issue.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.