EulerLawson fixed time step parameter

New to Julia
1

I’m trying to use the exponential runga kutta methods but I’m recieving the following error

LoadError: Fixed timestep methods require a choice of dt or choosing the tstops

When I go to supply the keyword dt = .1 to solve I get the error

LoadError: type ODEFunction has no field f1

Here’s the stack trace for the second error:

Stacktrace:
  [1] getproperty
    @ ./Base.jl:38 [inlined]
  [2] alg_cache(alg::LawsonEuler{10, true, Val{:forward}, true, nothing}, u::Matrix{Float64}, rate_prototype::Matrix{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float32}, uprev::Matrix{Float64}, uprev2::Matrix{Float64}, f::ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float32, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/twobl/src/caches/linear_nonlinear_caches.jl:194
  [3] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, alg::LawsonEuler{10, true, Val{:forward}, true, nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float32, force_dtmin::Bool, adaptive::Bool, gamma::Int64, abstol::Nothing, reltol::Nothing, qmin::Int64, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Int64, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/twobl/src/solve.jl:312
  [4] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/twobl/src/solve.jl:5 [inlined]
  [5] #solve_call#24
    @ ~/.julia/packages/DiffEqBase/WDtlB/src/solve.jl:451 [inlined]
  [6] solve_up(prob::ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, sensealg::Nothing, u0::Matrix{Float64}, p::SciMLBase.NullParameters, args::LawsonEuler{0, true, Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/WDtlB/src/solve.jl:814
  [7] #solve#29
    @ ~/.julia/packages/DiffEqBase/WDtlB/src/solve.jl:784 [inlined]
  [8] batch_func(i::Int64, prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::LawsonEuler{0, true, Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/mgXTt/src/ensemble/basic_ensemble_solve.jl:92
  [9] #469
    @ ~/.julia/packages/SciMLBase/mgXTt/src/ensemble/basic_ensemble_solve.jl:129 [inlined]
 [10] (::Base.var"#928#933"{SciMLBase.var"#469#470"{Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}}, EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, LawsonEuler{0, true, Val{:forward}, true, nothing}}})(r::Base.RefValue{Any}, args::Tuple{Int64})
    @ Base ./asyncmap.jl:100
 [11] macro expansion
    @ ./asyncmap.jl:234 [inlined]
 [12] (::Base.var"#944#945"{Base.var"#928#933"{SciMLBase.var"#469#470"{Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}}, EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, LawsonEuler{0, true, Val{:forward}, true, nothing}}}, Channel{Any}, Nothing})()
    @ Base ./task.jl:484
Stacktrace:
  [1] (::Base.var"#938#940")(x::Task)
    @ Base ./asyncmap.jl:177
  [2] foreach(f::Base.var"#938#940", itr::Vector{Any})
    @ Base ./abstractarray.jl:2774
  [3] maptwice(wrapped_f::Function, chnl::Channel{Any}, worker_tasks::Vector{Any}, c::UnitRange{Int64})
    @ Base ./asyncmap.jl:177
  [4] wrap_n_exec_twice
    @ ./asyncmap.jl:153 [inlined]
  [5] #async_usemap#923
    @ ./asyncmap.jl:103 [inlined]
  [6] #asyncmap#922
    @ ./asyncmap.jl:81 [inlined]
  [7] pmap(f::Function, p::Distributed.CachingPool, c::UnitRange{Int64}; distributed::Bool, batch_size::Int64, on_error::Nothing, retry_delays::Vector{Any}, retry_check::Nothing)
    @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/pmap.jl:126
  [8] solve_batch(prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::LawsonEuler{0, true, Val{:forward}, true, nothing}, ensemblealg::EnsembleDistributed, II::UnitRange{Int64}, pmap_batch_size::Int64; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/mgXTt/src/ensemble/basic_ensemble_solve.jl:128
  [9] macro expansion
    @ ./timing.jl:383 [inlined]
 [10] __solve(prob::EnsembleProblem{ODEProblem{Matrix{Float64}, Tuple{Float32, Float32}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, 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}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::LawsonEuler{0, true, Val{:forward}, true, nothing}, ensemblealg::EnsembleDistributed; trajectories::Int64, batch_size::Int64, pmap_batch_size::Int64, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :maxiters, :dt), Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Float64}}})
    @ SciMLBase ~/.julia/packages/SciMLBase/mgXTt/src/ensemble/basic_ensemble_solve.jl:56
 [11] #solve#31
    @ ~/.julia/packages/DiffEqBase/WDtlB/src/solve.jl:830 [inlined]
 [12] top-level scope
    @ ./timing.jl:263
 [13] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [14] top-level scope
    @ REPL[1]:1
in expression starting at /Users/marder/Dropbox/work/marder/homeo/runSimulations_LONG.jl:215

I also tried putting dt = .1 in ODEproblem but got the same error. Ideas about what’s happening here? thanks!

2

As documented, it’s an argument of solve.

https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/

3

I tried this, see line 3 of this post! I get another error.

4

What code did you run?

5

It’s a bit long but here it is:

using DifferentialEquations, Random, MAT

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

## Reversal Potentials
EK = -80.0;
EH = -20.0;
ENa = 50.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) = sigmd(V,0, 25.5, -5.29);  # m^3
NaHinf(V) = sigmd(V,0, 48.9, 5.18);  # h
CaTMinf(V) = sigmd(V,0, 27.1, -7.20);  # m^3
CaTHinf(V) = sigmd(V,0, 32.1, 5.50);  # h
CaSMinf(V) = sigmd(V,0, 33.0, -8.1);  # m^3
CaSHinf(V) = sigmd(V,0, 60.0, 6.20);  # h
HMinf(V) = sigmd(V,0, 70.0, 6.0);  # m
KdMinf(V) = sigmd(V,0, 12.3, -11.8);  # m^4
KCaMinf(V,IntCa) = (IntCa/(IntCa + 3.0))*sigmd(V,0, 28.3, -12.6); # m^4
AMinf(V) = sigmd(V,0, 27.2, -8.70); # m^3
AHinf(V) = sigmd(V,0, 56.9, 4.90);  # h

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

# Calcium Reverse Potential
R = 8.314*1000;  # Ideal Gas Constant (*10^3 to put into mV)
temp = 10;  # 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(m,h,V) = iIonic(7.85160170e1, m, h, 3.0, V, ENa);
iCaT(m,h,V,CaIn) = iIonic(3.36541720e1, m, h, 3.0, V, ECa(CaIn));
iCaS(m,h,V,CaIn) = iIonic(3.81677140e0, m, h, 3.0, V, ECa(CaIn));
iH(m,V) = iIonic(9.28628921e-4,  m, 1.0, 1.0, V, EH);
iKd(m,V) = iIonic(8.73559027e1, m, 1.0, 4.0, V, EK);
iKCa(m,V) = iIonic(4.76282766e1, m, 1.0, 4.0, V, EK);
iA(m,h,V) = iIonic(2.92172291e1, m, h, 3.0, V, EK);
iL(V) = iIonic(.01, 1.0, 1.0, 1.0, V, EL);
iApp(t) = 0;


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

function f(du,u,p,t)
    du[1] = -iL(u[1])-iNa(u[2],u[3],u[1])-iCaT(u[4],u[5],u[1],u[13])-iCaS(u[6],u[7],u[1],u[13])-iH(u[8],u[1])-iKd(u[9],u[1])-iKCa(u[10],u[1])-iA(u[11],u[12],u[1]) + iApp(t)
    du[2]  = (NaMinf(u[1]) - u[2])/tauNaM(u[1])
    du[3]  = (NaHinf(u[1]) - u[3])/tauNaH(u[1])
    du[4]  = (CaTMinf(u[1]) - u[4])/tauCaTM(u[1])
    du[5]  = (CaTHinf(u[1]) - u[5])/tauCaTH(u[1])
    du[6]  = (CaSMinf(u[1]) - u[6])/tauCaSM(u[1])
    du[7]  = (CaSHinf(u[1]) - u[7])/tauCaSH(u[1])
    du[8]  = (HMinf(u[1]) - u[8])/tauHM(u[1])
    du[9]  = (KdMinf(u[1]) - u[9])/tauKdM(u[1])
    du[10] = (KCaMinf(u[1],u[13]) - u[10])/tauKCaM(u[1])
    du[11] = (AMinf(u[1]) - u[11])/tauAM(u[1])
    du[12] = (AHinf(u[1]) - u[12])/tauAH(u[1])
    du[13] = (-.94*(iCaT(u[4],u[5],u[1],u[13])+iCaS(u[6],u[7],u[1],u[13]))-u[13]+.05)/20
end

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

ICs = [-60.0 rand(Float64, (1, 11)) 50*rand()];

tspan = (0.0f0,3000.0f0);
prob = ODEProblem(f,ICs,tspan);
function prob_func(prob,i,repeat)
  remake(prob,u0=ICs)
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sim = solve(ensemble_prob,LawsonEuler(),EnsembleDistributed(),trajectories=1, maxiters=1e12, dt=.1)

############################################################################################
### SAVING
## loop through trajectories lol

trj = 1

hash = randstring(['A':'Z'; '0':'9'], 6);
print(hash)
prfx = string(hash,"_AlonsoBnds_LiuAMat_TGT.mat");
arraySim = Array(sim);

# Trace
V_time = sim[trj].t[begin:5:end];
V = arraySim[1,1,begin:5:end,trj];

outFile = matopen(prfx,"w");
write(outFile,"V_time",V_time);
write(outFile,"V",V);
close(outFile);
6

Was someone able to replicate the error that I’m getting? Thanks!

7

It seems that LawsonEuler is a method for SplitODEProblems, and you are constructing an ODEProblem.

Basically, it is suited for systems of the form

\frac{du}{dt} = f_1(u, t) + f_2(u, t),

where typically one of those functions is represents a numerically stiff term. Are you sure you want to use that scheme to solve your problem?

1 Like
8

Yes that would be the error, the method is only for SplitODEProblems. (It could throw a better error message though, that’s something we should fix).

1 Like
9

okay ill take a look, thanks!