JuliaApproximation / ApproxFun.jl

Julia package for function approximation
http://juliaapproximation.github.io/ApproxFun.jl/
Other
538 stars 71 forks source link

Rounding errors lead to empty domain error when using abs #906

Open macd opened 1 year ago

macd commented 1 year ago

Because of rounding (I assume), when calculating roots we do not recover, say, the domain endpoint but something very close, ie 0.9999999999999994 instead of 1.0. If we try to take the abs of that Fun, it will throw an error complaining that the "Domain cannot be empty". The following shows both the bad behavior and, by perturbing one of the coefficients, the good behavior.

macd@macd-NUC9:~$ julia19
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.9.2 (2023-07-05)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia @v1.9> using ApproxFun

julia @v1.9> rcoefs = [7//32, -23//64, 9//32, -9//64]
4-element Vector{Rational{Int64}}:
   7//32
 -23//64
   9//32
  -9//64

julia @v1.9> f1 = Fun(Chebyshev(), Float64.(rcoefs))
Fun(Chebyshev(), [0.21875, -0.359375, 0.28125, -0.140625])

julia @v1.9> f2 = Fun(Chebyshev(), Float64.(rcoefs))
Fun(Chebyshev(), [0.21875, -0.359375, 0.28125, -0.140625])

julia @v1.9> f2.coefficients[2] += eps(1.0)
-0.3593749999999998

julia @v1.9> roots(f1)
3-element Vector{Float64}:
 -0.33333333333333326
  0.3333333333333333
  0.9999999999999994

julia @v1.9> roots(f2)
3-element Vector{Float64}:
 -0.3333333333333335
  0.33333333333333265
  1.0

julia @v1.9> abs(f2)
Fun(ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}}(PiecewiseSegment{Float64, Vector{Float64}}([-1.0, -0.3333333333333335, 0.33333333333333265, 1.0])), [1.0, 1.11022e-16, 8.67362e-17, 1.91687e-16, 0.09375, -0.03125, -0.03125, -0.00520833, 0.00520833, -0.00520833])

julia @v1.9> abs(f1)
ERROR: ArgumentError: Domain cannot be empty
Stacktrace:
  [1] Chebyshev
    @ ~/julia-versions/dot_julia/packages/ApproxFunOrthogonalPolynomials/hVSG5/src/Spaces/Chebyshev/Chebyshev.jl:16 [inlined]
  [2] Chebyshev
    @ ~/julia-versions/dot_julia/packages/ApproxFunOrthogonalPolynomials/hVSG5/src/Spaces/Chebyshev/Chebyshev.jl:22 [inlined]
  [3] canonicalspace
    @ ~/julia-versions/dot_julia/packages/ApproxFunOrthogonalPolynomials/hVSG5/src/Spaces/Ultraspherical/DirichletSpace.jl:30 [inlined]
  [4] checkcanonicalspace
    @ ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Space.jl:477 [inlined]
  [5] CanonicalTransformPlan
    @ ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Space.jl:484 [inlined]
  [6] plan_transform
    @ ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Space.jl:487 [inlined]
  [7] transform(S::ChebyshevDirichlet{1, 1, Segment{Float64}, Float64}, vals::Vector{Float64})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Space.jl:532
  [8] *(P::ApproxFunBase.TransformPlan{Float64, ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}}, false, Nothing}, vals::Vector{Float64})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Spaces/ContinuousSpace.jl:73
  [9] transform(S::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}}, vals::Vector{Float64})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/Space.jl:532
 [10] default_Fun(T::Type, f::Function, d::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}}, pts::Vector{Float64}, shouldsplat::Val{false})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:48
 [11] default_Fun(f::ComposedFunction{typeof(abs), Fun{Chebyshev{ChebyshevInterval{Float64}, Float64}, Float64, Vector{Float64}}}, d::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}}, n::Int64, shouldsplat::Val{false})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:62
 [12] _default_Fun(f::ComposedFunction{typeof(abs), Fun{Chebyshev{ChebyshevInterval{Float64}, Float64}, Float64, Vector{Float64}}}, d::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:117
 [13] default_Fun(f::Function, d::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:96
 [14] Fun(f::Function, d::ContinuousSpace{Float64, Float64, PiecewiseSegment{Float64, Vector{Float64}}})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:225
 [15] Fun(f::Function, d::PiecewiseSegment{Float64, Vector{Float64}})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/constructors.jl:205
 [16] splitmap(g::Function, d::ChebyshevInterval{Float64}, pts::Vector{Float64})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/specialfunctions.jl:2
 [17] abs(f::Fun{Chebyshev{ChebyshevInterval{Float64}, Float64}, Float64, Vector{Float64}})
    @ ApproxFunBase ~/julia-versions/dot_julia/packages/ApproxFunBase/rVDQs/src/specialfunctions.jl:31
 [18] top-level scope
    @ REPL[9]:1

julia @v1.9>