MethodError: no method matching similar(::Float64)

andreasvarga commented 5 months ago

The following code (also used in #57) fails:

using OrdinaryDiffEq
using IRKGaussLegendre

function x(t)
    u0 = [1.]
    tspan = (1.,t)
    prob = ODEProblem(Ric1ODE!, u0, tspan)
    sol = solve(prob, IRKGaussLegendre.IRKGL16(maxtrials=4); adaptive = true, reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)
    return sol(t)
end

function Ric1ODE!(du, u, pars, t) 
    du[1] = -2u[1]-1+u[1]^2
end

f(u, p, t) = x(t)[1]^2
u0 = 0.
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
# sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8); 
sol = solve(prob, IRKGaussLegendre.IRKGL16(maxtrials=4), reltol = 1e-8, abstol = 1e-8); 

sol(1)


ERROR: MethodError: no method matching similar(::Float64)

Closest candidates are:
  similar(::Union{LinearAlgebra.Adjoint{T, var"#s971"}, LinearAlgebra.Transpose{T, var"#s971"}} where {T, var"#s971"<:(AbstractVector)})
   @ LinearAlgebra C:\Users\Andreas\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\adjtrans.jl:329
  similar(::Union{LinearAlgebra.Adjoint{T, var"#s971"}, LinearAlgebra.Transpose{T, var"#s971"}} where {T, var"#s971"<:(AbstractVector)}, ::Type{T}) where T
   @ LinearAlgebra C:\Users\Andreas\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\adjtrans.jl:330
  similar(::Union{LinearAlgebra.Adjoint{T, S}, LinearAlgebra.Transpose{T, S}} where {T, S})
   @ LinearAlgebra C:\Users\Andreas\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\adjtrans.jl:333
  ...

Stacktrace:
 [1] __solve(::ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::IRKGL16{1, 4, true, false, false, false, Float64, 6}; dt::Float64, maxiters::Int64, save_everystep::Bool, adaptive::Bool, reltol::Float64, abstol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ IRKGaussLegendre C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16Solver.jl:190
 [2] __solve
   @ C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16Solver.jl:103 [inlined]
 [3] #solve_call#34
   @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:608 [inlined]
 [4] solve_call
   @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:566 [inlined]
 [5] #solve_up#42
   @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1057 [inlined]
 [6] solve_up
   @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1043 [inlined]
 [7] #solve#40
   @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:980 [inlined]
 [8] top-level scope
   @ c:\Users\Andreas\Documents\software\Julia\PeriodicSystems.jl\test\test_ODE.jl:26

The alternative call (to be commented out) to Tsit5() works correctly and the result is 3.3006156893224454.

ChrisRackauckas commented 5 months ago

This method doesn't handle scalar ODEs, i.e. u0 = 0., but u0 = [0.] should work.

andreasvarga commented 5 months ago

I redefined the functions in the outer integration as follows:

f(u, p, t) = x(t)[1]^2
function f1!(du, u, pars, t)
    du[1] = pars(t)
end
u0 = [0.]; u00 = 0.;
tspan = (0.0, 1.0)
pars = t->x(t)[1]^2
prob = ODEProblem(f, u00, tspan)
prob1 = ODEProblem(f1!, u0, tspan, pars)

and obtained the following results:

julia> sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8); sol(1)
3.3006156893224454

julia> sol1 = solve(prob1, Tsit5(), reltol = 1e-8, abstol = 1e-8); sol1(1)[1]
3.3006156893224454

julia> sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(maxtrials=4), reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
ERROR: NaN tspan is set in the problem or chosen in the init/solve call.
Note that -Inf and Inf values are allowed in the timespan for solves
which are terminated via callbacks, however NaN values are not allowed
since the direction of time is undetermined.

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Stacktrace:
  [1] get_concrete_tspan(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Ric1ODE!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool, kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:u0, :p, :adaptive, :reltol, :abstol, :save_everystep), Tuple{Vector{Float64}, SciMLBase.NullParameters, Bool, Float64, Float64, Bool}}}, p::SciMLBase.NullParameters)
    @ DiffEqBase C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1264
  [2] get_concrete_problem(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(Ric1ODE!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool; kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:u0, :p, :adaptive, :reltol, :abstol, :save_everystep), Tuple{Vector{Float64}, SciMLBase.NullParameters, Bool, Float64, Float64, Bool}}})
    @ DiffEqBase C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1146
  [3] get_concrete_problem
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1144 [inlined]
  [4] #solve_up#42
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1051 [inlined]
  [5] solve_up
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1043 [inlined]
  [6] #solve#40
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:980 [inlined]
  [7] x(t::Float64)
    @ Main c:\Users\Andreas\Documents\software\Julia\PeriodicSystems.jl\test\test_ODE.jl:8
  [8] (::var"#7#8")(t::Float64)
    @ Main c:\Users\Andreas\Documents\software\Julia\PeriodicSystems.jl\test\test_ODE.jl:27
  [9] f1!(du::Vector{Float64}, u::Vector{Float64}, pars::var"#7#8", t::Float64)
    @ Main c:\Users\Andreas\Documents\software\Julia\PeriodicSystems.jl\test\test_ODE.jl:23
 [10] (::SciMLBase.Void{typeof(f1!)})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\Andreas\.julia\packages\SciMLBase\8XHkk\src\utils.jl:481
 [11] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{typeof(f1!)}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::Function, arg4::Float64)
    @ FunctionWrappers C:\Users\Andreas\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
 [12] macro expansion
    @ C:\Users\Andreas\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
 [13] do_ccall
    @ C:\Users\Andreas\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
 [14] FunctionWrapper
    @ C:\Users\Andreas\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
 [15] _call
    @ C:\Users\Andreas\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
 [16] FunctionWrappersWrapper
    @ C:\Users\Andreas\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
 [17] ODEFunction
    @ C:\Users\Andreas\.julia\packages\SciMLBase\8XHkk\src\scimlfunctions.jl:2407 [inlined]
 [18] IRKstep_adaptive!(s::Int64, j::Int64, ttj::Vector{Float64}, tf::Float64, uj::Vector{Float64}, ej::Vector{Float64}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, var"#7#8", ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, var"#7#8", Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, var"#7#8", Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, var"#7#8", ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, var"#7#8", ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, dts::Vector{Float64}, coeffs::tcoeffs{Float64}, cache::IRKGaussLegendre.tcache{Vector{Float64}, Vector{Float64}, Float64, Vector{Float64}, Float64}, maxiters::Int64, maxtrials::Int64, initial_interp::Bool, abstol::Float64, reltol::Float64, adaptive::Bool, threading::Bool)
    @ IRKGaussLegendre C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16step_adaptive_seq.jl:92
 [19] IRKStep_seq!
    @ C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16Solver.jl:545 [inlined]
 [20] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, var"#7#8", ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, var"#7#8", Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, var"#7#8", Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, var"#7#8", ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, var"#7#8", ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::IRKGL16{1, 4, true, false, false, false, Float64, 6}; dt::Float64, maxiters::Int64, save_everystep::Bool, adaptive::Bool, reltol::Float64, abstol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ IRKGaussLegendre C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16Solver.jl:434
 [21] __solve
    @ C:\Users\Andreas\.julia\packages\IRKGaussLegendre\4pr2x\src\IRKGL16Solver.jl:103 [inlined]
 [22] #solve_call#34
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:608 [inlined]
 [23] solve_call
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:566 [inlined]
 [24] #solve_up#42
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1057 [inlined]
 [25] solve_up
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:1043 [inlined]
 [26] #solve#40
    @ C:\Users\Andreas\.julia\packages\DiffEqBase\7s9cb\src\solve.jl:980 [inlined]
 [27] top-level scope
    @ REPL[36]:1

Is this my error or is something wrong within IRKGaussLegendre ?

andreasvarga commented 5 months ago

With the function definition

function x(t)
    u0 = [1.]
    tspan = (1.,t)
    prob = ODEProblem(Ric1ODE!, u0, tspan)
    sol = solve(prob, Tsit5(); reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)
    return sol(t)
end

the following call works on Julia 1.9.2 (but not on Julia 1.10).

julia> sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
3.3006156872910055

mikelehu commented 5 months ago

Hi,

If you specify a initial step size, the following call works on both Julia 1.9.1 and Julia 1.10.0:

julia> sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt=1., reltol = 1e-8, abstol = 1e-8); sol2(1)[1] 3.3006156872882713

andreasvarga commented 5 months ago

Thanks for fixing that case. However, I still have problems understanding the interface logic of this function regarding the choice and role of keyword parameters like dt and adaptive. For example, on Julia 1.9.2:

sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.0000001,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]

does not work, but

sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.000001,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.0,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), adaptive=true, reltol = 1e-8, abstol = 1e-8); sol2(1)[1]

all work. On Julia 1.10, none of the following work

sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.0000001,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.000001,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 0.0,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), adaptive=true, reltol = 1e-8, abstol = 1e-8); sol2(1)[1]

but

sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt=0.00001, reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), adaptive=true, dt=1, reltol = 1e-8, abstol = 1e-8); sol2(1)[1]
sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt=1, reltol = 1e-8, abstol = 1e-8); sol2(1)[1]

all work.

Sorry for insisting on these issues, but this software is aimed to serve as the default option for solving periodic Riccati differential equations in my package PeriodicSystems, and I would like very much to master all provided facilities.

ChrisRackauckas commented 3 months ago

Solved in https://github.com/SciML/IRKGaussLegendre.jl/pull/64 by a modification that makes the error estimate of zero behave like machine epsilon.

mikelehu commented 3 months ago

I'm sorry for the delay, but I have been real busy. Now the code works correctly, I have tested it with all the cases that you indicated in the issue 61. I have added a new notebook with those tests (Examples/Riccati Differencial Equations.ipynb)
There was a bug in our implementation of the adaptive method (not related to the Julia version)
About "adaptive" and "dt" arguments:

1) “adaptive” argument can take two values: adaptive=true (default value), for adaptive timestepping. adaptive=false, for constant timestepping.

2) “dt” is the size of the initial step of numerical integration. In the case that the user does not specify one, the code makes this choice.

For constant step integrations, this dt will be used for all steps of the integration; in adaptive step integrations it will be adjusted according to the tolerance indicated by the user.

Regards, Mikel

SciML / IRKGaussLegendre.jl

MethodError: no method matching similar(::Float64) #61