Troubleshooting Integration

New to Julia
1

Hi,

I’m having trouble thinking about how to troubleshoot the following equation. I’ve isolated the time step that produces the error and have examined many values that the error suggests to look at (for instance, the values of ECa, Ca, ICaT, ICaS, dCa/dt - equation 13, etc). Nothing looks a miss. The next step for me was to examine the values the integrator is producing. But I’m having trouble using the integrator interface. Below is the code and things I’ve tried. Do people see what’s going on? Thanks!

Code

using DifferentialEquations

endd = 99550.9  # last timestep that works

#endd = 99551.0 # first timestep that fails

function EK(t)
    if t < 0
        return -80
    elseif t >= 0 && t < 50000.0
        return -80
    elseif t >= 50000.0 && t < 150000.0
        return -80
    else
        return -80
    end
end

#####################################################################
### Parameters

## Sensor Parameters
alphaOffset = 0.075 # alpha in between
GF =  53;
GS =  3;
GD =  1.0;
gamma = 10e-7;
delta = 10e-4;
sensorFuseDelta = .005;

## Targets
FBar = .25;
SBar = .03; 
DBar = .02;

## Homeostatic Time Scales
tauG = 4000.0;
tauHalfac = Inf;
tauS = 2000.0;
tauAlpha = 2000.0;

## Reversal Potentials
EH = -20.0;
ENa = 30.0;
EL = -50.0;

#####################################################################

### Define generic functions
## Generic Activation/Inactication Curves
sigmd(x,translate,center,reactance) = 1/(1+exp(((x-translate) + center)/reactance));
## Define the time constant function for ion channels
tauX(Volt, CT, DT, AT, BT, halfAc) = CT - DT/(1 + exp(((Volt-halfAc) + AT)/BT));
## Some intrinsic currents require a special time constant function
spectau(Volt, CT, DT, AT, BT, AT2, BT2, halfAc) = CT + DT/(exp(((Volt-halfAc) + AT)/BT) + exp(((Volt-halfAc) + AT2)/BT2));
## Define the ionic currents; q is the exponent of the activation variable m
iIonic(g, m, h, q, Volt, Erev) = g*(m^q)*h*(Volt - Erev);

#####################################################################

# Ion Channel Acitvation/Inactivation Curves
NaMinf(V,trns) = sigmd(V,trns, 25.5, -5.29);  # m^3
NaHinf(V,trns) = sigmd(V,trns, 48.9, 5.18);  # h
CaTMinf(V,trns) = sigmd(V,trns, 27.1, -7.20);  # m^3
CaTHinf(V,trns) = sigmd(V,trns, 32.1, 5.50);  # h
CaSMinf(V,trns) = sigmd(V,trns, 33.0, -8.1);  # m^3
CaSHinf(V,trns) = sigmd(V,trns, 60.0, 6.20);  # h
HMinf(V,trns) = sigmd(V,trns, 70.0, 6.0);  # m
KdMinf(V,trns) = sigmd(V,trns, 12.3, -11.8);  # m^4
KCaMinf(V,trns,IntCa) = (IntCa/(IntCa + 3.0))*sigmd(V,trns, 28.3, -12.6); # m^4
AMinf(V,trns) = sigmd(V,trns, 27.2, -8.70); # m^3
AHinf(V,trns) = sigmd(V,trns, 56.9, 4.90);  # h

# Time Constants (ms)
tauNaM(V,trns) = tauX(V, 1.32, 1.26, 120.0, -25.0, trns);
tauNaH(V,trns) = tauX(V, 0.00, -0.67, 62.9, -10.0, trns)*tauX(V, 1.50, -1.00, 34.9, 3.60, trns);
tauCaTM(V,trns) = tauX(V, 21.7, 21.3, 68.1, -20.5, trns);
tauCaTH(V,trns) = tauX(V, 105.0, 89.8, 55.0, -16.9, trns);
tauCaSM(V,trns) = spectau(V, 1.40, 7.00, 27.0, 10.0, 70.0, -13.0, trns);
tauCaSH(V,trns) = spectau(V, 60.0, 150.0, 55.0, 9.00, 65.0, -16.0, trns);
tauHM(V,trns) = tauX(V, 272.0, -1499.0, 42.2, -8.73, trns);
tauKdM(V,trns) = tauX(V, 7.20, 6.40, 28.3, -19.2, trns);
tauKCaM(V,trns) = tauX(V, 90.3, 75.1, 46.0, -22.7, trns);
tauAM(V,trns) = tauX(V, 11.6, 10.4, 32.9, -15.2, trns);
tauAH(V,trns) = tauX(V, 38.6, 29.2, 38.9, -26.5, trns);

# Calcium Reverse Potential
R = 8.314*1000.0;  # Ideal Gas Constant (*10^3 to put into mV)
temp = 10.0;  # Temperature in Celcius; Temperature of the Sea lawl
z = 2.0;  # Valence of Calcium Ions
Far = 96485.33;  # Faraday's Constant
CaOut = 3000.0;  # Outer Ca Concentration (uM)
ECa(CaIn) = ((R*(273.15 + temp))/(z*Far))*log(CaOut/CaIn);

# Ionic Currents (mV / ms)
iNa(g,m,h,V) = iIonic(g, m, h, 3.0, V, ENa);
iCaT(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iCaS(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iH(g,m,V) = iIonic(g,  m, 1.0, 1.0, V, EH);
iKd(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iKCa(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iA(g,m,h,V,t) = iIonic(g, m, h, 3.0, V, EK(t));
iL(V) = iIonic(.01, 1.0, 1.0, 1.0, V, EL);
iApp(t) = 0;

# Sensor Equations
F(FM,FH) = GF*FM^2*FH;
S(SM,SH) = GS*SM^2*SH;
D(DM)    = GD*DM^2; 
sensorFuse(EF,ES,ED) = exp(-abs(EF)/(sensorFuseDelta*10))*exp(-abs(ES)/sensorFuseDelta)*exp(-abs(ED)/sensorFuseDelta);

############################################################################################

function f(du,u,p,t)
    du[1] = -iL(u[1])-iNa(u[14],u[2],u[3],u[1])-iCaT(u[15],u[4],u[5],u[1],u[13])-iCaS(u[16],u[6],u[7],u[1],u[13])-iH(u[17],u[8],u[1])-iKd(u[18],u[9],u[1],t)-iKCa(u[19],u[10],u[1],t)-iA(u[20],u[11],u[12],u[1],t)
    du[2]  = (NaMinf(u[1],u[21]) - u[2])/tauNaM(u[1],u[21])
    du[3]  = (NaHinf(u[1],u[22]) - u[3])/tauNaH(u[1],u[22])
    du[4]  = (CaTMinf(u[1],u[23]) - u[4])/tauCaTM(u[1],u[23])
    du[5]  = (CaTHinf(u[1],u[24]) - u[5])/tauCaTH(u[1],u[24])
    du[6]  = (CaSMinf(u[1],u[25]) - u[6])/tauCaSM(u[1],u[25])
    du[7]  = (CaSHinf(u[1],u[26]) - u[7])/tauCaSH(u[1],u[26])
    du[8]  = (HMinf(u[1],u[27]) - u[8])/tauHM(u[1],u[27])
    du[9]  = (KdMinf(u[1],u[28]) - u[9])/tauKdM(u[1],u[28])
    du[10] = (KCaMinf(u[1],u[29],u[13]) - u[10])/tauKCaM(u[1],u[29])
    du[11] = (AMinf(u[1],u[30]) - u[11])/tauAM(u[1],u[30])
    du[12] = (AHinf(u[1],u[31]) - u[12])/tauAH(u[1],u[31])
    du[13] = (-.94*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]))-u[13]+.05)/20
    du[14] = ((1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[14]-gamma*u[14]^3)*u[40]/tauG
    du[15] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[15]-gamma*u[15]^3)*u[40]/tauG
    du[16] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[16]-gamma*u[16]^3)*u[40]/tauG
    du[17] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))*u[17]-gamma*u[17]^3)*u[40]/tauG
    du[18] = ((1*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[18]-gamma*u[18]^3)*u[40]/tauG
    du[19] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))*u[19]-gamma*u[19]^3)*u[40]/tauG
    du[20] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))*u[20]-gamma*u[20]^3)*u[40]/tauG
    du[21] = ((-1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[21]^3)*u[40]/tauHalfac
    du[22] = ((1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[22]^3)*u[40]/tauHalfac
    du[23] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[23]^3)*u[40]/tauHalfac
    du[24] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[24]^3)*u[40]/tauHalfac
    du[25] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[25]^3)*u[40]/tauHalfac
    du[26] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[26]^3)*u[40]/tauHalfac
    du[27] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))-delta*u[27]^3)*u[40]/tauHalfac
    du[28] = ((-1*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[28]^3)*u[40]/tauHalfac
    du[29] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))-delta*u[29]^3)*u[40]/tauHalfac
    du[30] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))-delta*u[30]^3)*u[40]/tauHalfac
    du[31] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))-delta*u[31]^3)*u[40]/tauHalfac
    du[32] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,14.8,1) - u[32])/.5
    du[33] = (sigmd(-1*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13])),0,-9.8,1) - u[33])/1.5
    du[34] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,7.2,1) - u[34])/50
    du[35] = (sigmd(-1*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13])),0,-2.8,1) - u[35])/60
    du[36] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,3,1) - u[36])/500
    du[37] = ((FBar - F(u[32],u[33])) - u[37])/tauS
    du[38] = ((SBar - S(u[34],u[35])) - u[38])/tauS
    du[39] = ((DBar - D(u[36])) - u[39])/tauS
    du[40] = (sigmd(-1*sensorFuse(u[37],u[38],u[39]),0,alphaOffset,-.01)-u[40])/tauAlpha
end

############################################################################################

ICs = [-50.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.8320755796233367 0.046612330545369285 0.4039904957344005 0.7335731823371116 0.516374395070418 0.7517868565780351 0.9012425121331321 0.3795887585487998 -0.49230754975054447 -0.650925587096205 0.3021172776981724 -0.2266327552455598 -0.1602097765267385 0.2278426433885099 -0.11586687638108684 -0.5353572538829641 0.46427761017041913 -0.626348202330232 0.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0 0.5];

tspan = (0.0,endd);
SaveTimes = (0.0:.1:endd); #in milliseconds

prob = ODEProblem(f,ICs,tspan);
@time sim = solve(prob,Tsit5(), maxiters=1e10, saveat=SaveTimes);

When run for t = 99551.0 I obtain the following error:

ERROR: DomainError with -867.2167835365867:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:292
  [3] log(x::Float64)
    @ Base.Math ./special/log.jl:257
  [4] ECa(CaIn::Float64)
    @ Main ./REPL[123]:1
  [5] iCaT(g::Float64, m::Float64, h::Float64, V::Float64, CaIn::Float64)
    @ Main ./REPL[125]:1
  [6] f(du::Matrix{Float64}, u::Matrix{Float64}, p::SciMLBase.NullParameters, t::Float64)
    @ Main ./REPL[137]:2
  [7] ODEFunction
    @ ~/.julia/packages/SciMLBase/BC0Kw/src/scimlfunctions.jl:1622 [inlined]
  [8] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/perform_step/low_order_rk_perform_step.jl:709
  [9] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/perform_step/low_order_rk_perform_step.jl:675 [inlined]
 [10] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:477
 [11] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:5 [inlined]
 [12] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:437 [inlined]
 [13] solve_up(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Matrix{Float64}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:maxiters, :saveat), Tuple{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:780
 [14] #solve#33
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:760 [inlined]
 [15] top-level scope
    @ ./timing.jl:263 [inlined]
 [16] top-level scope
    @ ./REPL[142]:0

To use the integrator interface I add the following to the end of my code:

integrator = init(prob,Tsit5());
step!(integrator, 99550.9, stop_at_tdt = false)

And I receive the following error:

ERROR: MethodError: no method matching step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, ::Float64; stop_at_tdt=false)
Closest candidates are:
  step!(::SciMLBase.DEIntegrator, ::Any) at ~/.julia/packages/SciMLBase/BC0Kw/src/integrator_interface.jl:735 got unsupported keyword argument "stop_at_tdt"
  step!(::SciMLBase.DEIntegrator, ::Any, ::Any) at ~/.julia/packages/SciMLBase/BC0Kw/src/integrator_interface.jl:735 got unsupported keyword argument "stop_at_tdt"
  step!(::OrdinaryDiffEq.ODEIntegrator) at ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/iterator_interface.jl:1 got unsupported keyword argument "stop_at_tdt"
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[73]:1
2

That’s positional, not keyword.

step!(integrator, 99550.9, false)

3

I think the timestepping just gets too aggressive and takes such a big step that the argument to the log becomes negative. One nice thing to trying in this case is to replace Base.log with NaNMath.log, returning NaN instead of raising an exception allows the integrator to notice something has gone wrong and try a smaller timestep. This lets your simulation complete when I tried it. In addition, AutoTsit5(Rosenbrock23()) was substantially faster than Tsit5()

4

Thanks!

  1. I’m unsure how the integrator works – I want to return something at a resolution of “.1” (As I specify in saveat). Does the integrator use an adaptive time step then “interpolate” to this resolution?
  2. What does AutoTsit5(Rosenbrock23()) do?
5

Thanks!

6

Does the integrator use an adaptive time step then “interpolate” to this resolution?

Specifying saveat ensures that the solution has a point at this time, the solver is free to actually step here or interpolate as needed to achieve the tolerances in the most efficient way possible.

AutoTsit5(Rosenbrock23()) is a Composite Algorithm that chooses among several algorithms based on the detected stiffness of the problem.

1 Like