ApproxFun.jl icon indicating copy to clipboard operation
ApproxFun.jl copied to clipboard

Solving ODE system fails with MethodError depending on interval endpoint

Open vram0gh opened this issue 2 years ago • 1 comments

The attached code shows two cases of attempting to solve an ODE problem using ApproxFun. Changing the interval endpoint, without any other code changes, results in a shift from successful solutions to a MethodError. The code is a self-contained reproduction case that has been stripped down, with hardcoded initial values. It results in the error whose stacktrace is given below on an AMD64 system running Linux; the ApproxFun version used was 0.12.5.

Invoke works() to see a successful solve, and fail_methoderror() to see the error/stacktrace corresponding to:

ERROR: MethodError: no method matching 
setdomain(::PiecewiseSpace{Tuple{Chebyshev{ChebyshevInterval{Float64},Float64},JacobiWeight{Chebyshev{Interval{:closed,:closed,Float64},Float64},Interval{:closed,:closed,Float64},Float64,Int64},JacobiWeight{Chebyshev{Interval{:closed,:closed,Float64},Float64},Interval{:closed,:closed,Float64},Float64,Int64}},DomainSets.UnionDomain{Float64,Tuple{ChebyshevInterval{Float64},Interval{:closed,:closed,Float64},Interval{:closed,:closed,Float64}}},Float64}, ::Interval{:closed,:closed,Float64})
Closest candidates are:
  setdomain(::ApproxFunBase.SumSpace, ::Domain) at /home/user/.julia/packages/ApproxFunBase/Gra7W/src/Spaces/SumSpace.jl:105
  setdomain(::PiecewiseSpace, ::DomainSets.UnionDomain) at /home/user/.julia/packages/ApproxFunBase/Gra7W/src/Spaces/SumSpace.jl:108
  setdomain(::ApproxFunBase.ArraySpace, ::Domain) at /home/user/.julia/packages/ApproxFunBase/Gra7W/src/Spaces/ArraySpace.jl:67
  ...
Stacktrace:
 [1] setdomain(::Fun{PiecewiseSpace{Tuple{Chebyshev{ChebyshevInterval{Float64},Float64},JacobiWeight{Chebyshev{Interval{:closed,:closed,Float64},Float64},Interval{:closed,:closed,Float64},Float64,Int64},JacobiWeight{Chebyshev{Interval{:closed,:closed,Float64},Float64},Interval{:closed,:closed,Float64},Float64,Int64}},DomainSets.UnionDomain{Float64,Tuple{ChebyshevInterval{Float64},Interval{:closed,:closed,Float64},Interval{:closed,:closed,Float64}}},Float64},Float64,Array{Float64,1}}, ::Interval{:closed,:closed,Float64}) at /home/user/.julia/packages/ApproxFunBase/Gra7W/src/Fun.jl:184
 [2] /(::Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}}, ::Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}}) at /home/user/.julia/packages/ApproxFunBase/Gra7W/src/specialfunctions.jl:98
 [3] /(::ApproxFun.DualFun{Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}},PlusOperator{Float64,Tuple{Int64,Int64}}}, ::ApproxFun.DualFun{Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}},ApproxFunBase.ConstantTimesOperator{PlusOperator{Float64,Tuple{Int64,Int64}},Float64}}) at /home/user/.julia/packages/ApproxFun/GOMNa/src/Extras/autodifferentiation.jl:55
 [4] (::var"#N#1"{Float64,Float64,Float64,Float64})(::ApproxFun.DualFun{Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}},ConstantOperator{Float64,Chebyshev{Interval{:closed,:closed,Float64},Float64}}}, ::Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}}) at /home/user/research_code/approxfun_issue_cf_min.jl:14
 [5] newton(::Function, ::Array{Fun{Chebyshev{Interval{:closed,:closed,Float64},Float64},Float64,Array{Float64,1}},1}; maxiterations::Int64, tolerance::Float64) at /home/user/.julia/packages/ApproxFun/GOMNa/src/Extras/autodifferentiation.jl:137
 [6] approxfun_final_ode_step(::Float64) at /home/user/research_code/approxfun_issue_cf_min.jl:21
 [7] solve_test at /home/user/research_code/approxfun_issue_cf_min.jl:28 [inlined]
 [8] fail_methoderror() at /home/user/research_code/approxfun_issue_cf_min.jl:36
 [9] top-level scope at REPL[3]:1

Notes:

  1. I used a fairly coarse tolerance in the Newton solver, and still get this error.
  2. It takes a number of minutes for the failing case to arrive at the error. A lot of work of some kind is occurring.
  3. A singularity occurs in the ODE solution within the interval, in both the works() and fail_methoderror() cases; in fact, this is the motivation for using ApproxFun.jl rather than DifferentialEquations.jl. I only care about the solution before it reaches the singularity, but do not know a priori where the singularity is.
  4. If this happens to be a pathological set of inputs, it would be very helpful to get some kind of deterministic error/exception, or to be able to set limits on how hard ApproxFun tries to build suitable approximations (whether memory/size of matrices, tolerances, etc.). Since in an arbitrary problem solve, the location of the singularity is not known a priori , robustness would be needed in the solution attempt.

Code (I wasn't able to upload without appending a .txt extension): approxfun_issue_cf_min.jl.txt

vram0gh avatar Aug 22 '21 03:08 vram0gh

When dividing by a Fun, ApproxFun first transforms the domain to the canoncial domain of the associated space. If the divisor has roots, the domain is split at the roots and the space converted into a PiecewiseSpace. After the division, the domain of the space should be transformed back into to the original domain. This is where the issue occurs. It seems there currently is no implementation of setdomain which will remap the fragmented domain to the coordinates of the initial domain. As a temporary fix you can reparametrize your problem to live on the canonical domain. That fixed the problem for me.

LucasAschenbach avatar Aug 28 '23 20:08 LucasAschenbach