MethodError: no method matching similar(::Float64)

[andreasvarga]

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]

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

[andreasvarga]

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]

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]

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]

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]

Solved in #64 by a modification that makes the error estimate of zero behave like machine epsilon.

[mikelehu]

• 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