SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.82k stars 222 forks source link

Integration time/time-step sanity check #419

Closed daviehh closed 5 years ago

daviehh commented 5 years ago

ODE: around 500 coupled odes (written in complex-valued matrices), involves some random numbers. They are hierarchical equations, where higher orders are supposed to be close to zero. (unable to share full equation set for now, so this is just a general query)

[it is understood that the solver provided may fail for some problems; the main purpose of this issue is to discuss/report that some errors makes more logical sense for the ode solver to handle, rather than a NaN error reported by julia]

Solver: AutoVern8(Rosenbrock23()) or AutoTsit5(Rosenbrock23()), with autodiff=false (for complex numbers)

Problem encountered:

for problem function f(du, u,p,t), I got an error message that I'm accessing t=NaN. Turns out both t and u have become NaN, while du sometimes is NaN, sometimes is a very, very big number (20-digit or so?) while in reality it should be around 0~5.

on the other hand, when I manually convert the equation by splitting the complex into real and imaginary part, I normally get

┌ Warning: Instability detected. Aborting └ @ DiffEqBase ~/.julia/packages/DiffEqBase/VM2Bp/src/integrator_interface.jl:162

while occasionally also getting the accessing t at NaN error message from julia, full message at the end. The warning message and sol.retcode being not :Success is the "correct" behavior, while the t=NaN error comes from julia and seems to be not handle by the differential equation solver. Is there some checks in the solver to make sure that the timestep/final time does not go over tspan, or, handle NaN returned by the problem function (and return the error/warning as expected)? Thanks!

full error message:

ERROR: LoadError: BoundsError: attempt to access 601-element scale(interpolate(OffsetArray(::Array{Complex{Float64},1}, 0:602), BSpline(Cubic(Line(Interpolations.OnGrid())))), (0.000:0.050:30.000,)) with element type Complex{Float64} at index [NaN] Stacktrace: [1] throw_boundserror(::Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Tuple{Float64}) at ./abstractarray.jl:484 [2] ScaledInterpolation at /Users/.../.julia/packages/Interpolations/0OnYb/src/scaling/scaling.jl:66 [inlined] [3] dot_func(::Array{Float64,1}, ::Array{Float64,1}, ::NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}, ::Float64) at /Users/.../Downloads/tmp/diff/eqn_hfd_solver.jl:63 [4] ODEFunction at /Users/.../.julia/packages/DiffEqBase/VM2Bp/src/diffeqfunction.jl:106 [inlined] [5] perform_step!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},Array{Float64,1},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern8Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern8ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqDiffTools.UJacobianWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqBase.LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern8Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern8ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqDiffTools.UJacobianWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqBase.LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),DiffEqBase.CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}, ::OrdinaryDiffEq.Vern8Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern8ConstantCache{Float64,Float64}}, ::Bool) at /Users/.../.julia/packages/OrdinaryDiffEq/suHr2/src/perform_step/verner_rk_perform_step.jl:528 [6] perform_step! at /Users/.../.julia/packages/OrdinaryDiffEq/suHr2/src/perform_step/composite_perform_step.jl:43 [inlined] (repeats 2 times) [7] solve!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},Array{Float64,1},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern8Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern8ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqDiffTools.UJacobianWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqBase.LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern8Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern8ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqDiffTools.UJacobianWrapper{DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}}},DiffEqBase.LinSolveFactorize{typeof(lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),DiffEqBase.CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}) at /Users/.../.julia/packages/OrdinaryDiffEq/suHr2/src/solve.jl:329 [8] #__solve#457(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::DiffEqBase.ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/.../.julia/packages/OrdinaryDiffEq/suHr2/src/solve.jl:7 [9] __solve at /Users/.../.julia/packages/OrdinaryDiffEq/suHr2/src/solve.jl:6 [inlined] (repeats 3 times) [10] #solve#442(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::DiffEqBase.ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}) at /Users/.../.julia/packages/DiffEqBase/VM2Bp/src/solve.jl:39 [11] solve(::DiffEqBase.ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:zf, :zdt, :h, :l, :hdim, :pn, :a0, :gamma_eff),Tuple{Interpolations.ScaledInterpolation{Complex{Float64},1,Interpolations.BSplineInterpolation{Complex{Float64},1,OffsetArrays.OffsetArray{Complex{Float64},1,Array{Complex{Float64},1}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{Base.OneTo{Int64}}},Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},Float64,getfield(eqn_hfd_solver, Symbol("#h#152")){NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}},Array{Int64,2},Int64,Int64,Float64,Float64}},DiffEqBase.ODEFunction{true,typeof(eqn_hfd_solver.dot_func),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType}},OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Vern8,OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{typeof(lu!)},DataType},Rational{Int64},Float64}}) at /Users/.../.julia/packages/DiffEqBase/VM2Bp/src/solve.jl:27 [12] eqn_solve(::NamedTuple{(:pn, :tf, :zdt, :Gamma, :gamma_eff, :hsys, :l, :psi0, :frsn, :zn),Tuple{Int64,Float64,Float64,Int64,Float64,typeof(h),Array{Int64,2},Array{Float64,1},Int64,Int64}}) at /Users/.../Downloads/tmp/diff/eqn_hfd_solver.jl:185 [13] top-level scope at none:0

daviehh commented 5 years ago

also, looks like NaN returned by the problem function is handled correctly, e.g. with

function f(du, u,p,t)
    du = t>1 ? NaN : 0.1
end
ChrisRackauckas commented 5 years ago

It will be almost impossible to help you here without seeing any code. That said,

function f(du, u,p,t)
    du = t>1 ? NaN : 0.1
end

Note that the code you wrote there doesn't do anything. du = ... overrides the input, so no mutation occurs. Instead, you should be doing something like du[1] = ... or du .= ....

But a bigger question.

ODE: around 500 coupled odes (written in complex-valued matrices), involves some random numbers. They are hierarchical equations, where higher orders are supposed to be close to zero. (unable to share full equation set for now, so this is just a general query)

Are you sure it's an ODE? If you have random numbers in your derivative function, it is not an ODE. It might be a stochastic differential equation or a random ordinary differential equation. If you are taking new random numbers each time f is called, then the adaptive algorithm will sense that f is not defined and will reject the step, and keep doing so until the simulation ends or it goes unstable. Problems with stochasticity should use the SDE or RODE solvers for this purpose, or use adaptive=false to turn off adaptivity. But note that these methods are not guaranteed the same order of convergence when there's randomness anyways, whereas the SDE and RODE methods are made to converge (with higher order) in cases with randomness.

So anyways your issue is either because du isn't actually updated since you're shadowing instead of mutating it, or it's because of taking random numbers in your derivative function.

daviehh commented 5 years ago

Thanks! Sorry, that du = is a typo, in my actual code I have things like du[1:n1] = ...; the random number is generated and saved in p, then passes on to f (something like ODEProblem(f,u0,tspan,p)) so f is not changed every time it's called (no explicit random numbers in f). I'll look into using the SODE solvers, we are using a non-standard complex-valued noise process but it looks like there's a NoiseGrid function in the package that allows one to use arbitrarily noise process. Thanks! closing the issue now.