# Callbacks for parameters in a stiff system of ODEs?

**URL:** <https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580>\
**Category:** Modelling & Simulations\
**Created:** [June 17, 2020, 7:55am UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580 "2020-06-17T07:55:57Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![dynova](https://avatars.discourse-cdn.com/v4/letter/d/e47c2d/32.png) [@dynova](https://discourse.julialang.org/u/dynova)\
**Post date:** [June 17, 2020, 7:55am UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/1 "2020-06-17T07:55:57Z")

</div>

I have a stiff system of 21 ODEs, for which I want to use callback functions on some time-dependent parameters.

Looking for some help on this, as I’m running into some error messages and I’m new to callbacks. Not entirely sure how to incorporate the piecewise initial conditions.

Basically, I’m looking for:

```julia
v1 = 0 and v2 = 0 if t < a + d
v1 = gh and v2 = g(1-h) if t >= a+d
ydot = zeros(21,1) if t < a
ydot = ydot if t >= a

```

This is my code:

```julia
using LinearAlgebra, DifferentialEquations, MATLABDiffEq

a = 25;
d = 5;
T = 2323;
V = 1;
c = 0.1;
d = 0.35;
e = 0.5;
f = 0.2;
g = 1.0/3.0;
h = 0.3;
k = 0.12;
l = 1.0/11;
m = 0.6;
n = 0.1;
p = n*l;
q = (1-n)*l;
z = 0.01;
s = 0.2;
t = q;
u = 0.5*l;
v = 0.592;
w = v*u;
x = (1-v)*u;
y = 0.23245;
z = 13;

function RHS(ydot,y,par,t)

ydot = 0.0*zeros(21,1);
y = 0.0*zeros(21,1);
par = 0.0*zeros(23,1);
v1 = 0.0;
v2 = 0.0;
a = par[1];
d = par[2];
T = par[3];
d = par[4];
c = par[5];
e = par[6];
f = par[7];
g = par[8];
h = par[9];
k = par[10];
l = par[11];
m = par[12];
n = par[13];
p = par[14];
q = par[15];
z = par[16];
s = par[17];
t = par[18];
u = par[19];
w = par[20];
x = par[21];
y = par[22];
z = par[23];

fM = d * y[1] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
h = c * d * y[2] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;

ydot[1] = -fM - (v1*y[1]-v2*y[2]);
ydot[2] = -h + (v1*y[1]-v2*y[2]);
ydot[3] = fM - k*y[3]-(v1*y[3]-v2*y[4]);
ydot[4] = h - k*y[4]+(v1*y[3]-v2*y[4]);
ydot[5] = k*(y[3]-y[5])-z*y[5]-(v1*y[5]-v2*y[6]);
ydot[6] = k*(y[4]-y[6])-z*y[6]+(v1*y[5]-v2*y[6]);
ydot[7] = z*(y[5]+y[6])-k*y[7];
ydot[8] = k*(y[5]-y[8])-z*y[8]-(v1*y[8]-v2*y[9]);
ydot[9] = k*(y[6]-y[9])-z*y[9]+(v1*y[8]-v2*y[9]);
ydot[10] = z*(y[8]+y[9])-k*y[10] + k*y[7];
ydot[11] = k*(y[8]-y[11])-z*y[11]-(v1*y[11]-v2*y[12]);
ydot[12] = k*(y[9]-y[12])-z*y[12]+(v1*y[11]-v2*y[12]);
ydot[13] = z*(y[11]+y[12])-k*y[13] + k*y[10];
ydot[14] = k*(y[11]-y[14])-z*y[14]-(v1*y[14]-v2*y[15]);
ydot[15] = k*(y[12]-y[15])-z*y[15]+(v1*y[14]-v2*y[15]);
ydot[16] = z*(y[14]+y[15])-k*y[16] + k*y[13];
ydot[17] = m*k*y[14]-z*y[17]-(v1*y[17]-v2*y[18])-t*y[17];
ydot[18] = m*k*y[15]-z*y[18]+(v1*y[17]-v2*y[18])-t*y[18];
ydot[19] = (1-m)*k*y[14]-(z+s)*y[19]-(v1*y[19]-v2*y[20])-(p+q)*y[19];
ydot[20] = (1-m)*k*y[15]-(z+s)*y[20]+(v1*y[19]-v2*y[20])-(p+q)*y[20];
ydot[21] = (1-m)*k*(y[14]+y[15]);

return ydot
end

function condition(out,y,t,integrator) # Event when event_f(u,t) == 0
out[1] = integrator.t - integrator.par[1] - integrator.par[2];
out[2] = integrator.par[1] - integrator.t;
end

function affect!(integrator, idx)
if idx == 1
    v1 = integrator.par[8]*integrator.par[9];
    v2 = integrator.par[8]*(1-integrator.par[9]);
else
    ydot = 0.0*zeros(21,1); 
end
end

cb = VectorContinuousCallback(condition,affect!,2);
y0 = [2323 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1];
tspan = (0.0, 100.0);
prob = ODEProblem(RHS,y0,tspan);
sol = solve(prob, MATLABDiffEq.ode15s(), callback = cb, saveat = 1.0)

```

and the error with stack trace:

```julia
ERROR: LoadError: BoundsError
Stacktrace:
 [1] getindex at ./number.jl:78 [inlined]
 [2] RHS(::Array{Any,2}, ::Array{ModelingToolkit.Operation,2}, ::Array{Any,1}, ::ModelingToolkit.Operation) at /Users/dc/desktop/is/post.jl:60
 [3] (::ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Any,2}, ::Vararg{Any,N} where N) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/diffeqfunction.jl:248
 [4] modelingtoolkitize(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /Users/dc/.julia/packages/ModelingToolkit/wy693/src/systems/diffeqs/modelingtoolkitize.jl:19
 [5] __solve(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}; saveat::Float64, timeseries_errors::Bool, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}}}}) at /Users/dc/.julia/packages/MATLABDiffEq/b8nHE/src/MATLABDiffEq.jl:47
 [6] solve_call(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:callback, :saveat),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Float64}}}) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:96
 [7] #solve_up#454 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:122 [inlined]
 [8] #solve#453 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:106 [inlined]
 [9] top-level scope at /Users/dc/desktop/is/post.jl:106
 [10] include(::Module, ::String) at ./Base.jl:377
 [11] exec_options(::Base.JLOptions) at ./client.jl:288
 [12] _start() at ./client.jl:484
in expression starting at /Users/dc/desktop/is/post.jl:106

```

Thanks for any help!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 17, 2020, 11:01am UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/2 "2020-06-17T11:01:43Z")

</div>

There are a number of issues with your model. First off, MATLABDiffEq doesn’t allow for callbacks: it’ll throw an error when you get there anyways, so just use the normal methods (we cannot inject enough code into MATLAB to do it).

But secondly, the way you were referencing parameters wasn’t quite right. I think it might help if you use a LabelledArray, since then you can reference things by name. So I got a working code for you here:

```julia
using DifferentialEquations, LinearAlgebra, LabelledArrays

par = LVector(
    a = 25,
    d1 = 5,
    T = 2323,
    V = 1,
    c = 0.1,
    d2 = 0.35,
    e = 0.5,
    f = 0.2,
    g = 1.0/3.0,
    h = 0.3,
    k = 0.12,
    l = 1.0/11,
    m = 0.6,
    n = 0.1,
    p = n*l,
    q = (1-n)*l,
    z1 = 0.01,
    s = 0.2,
    t = q,
    u = 0.5*l,
    v = 0.592,
    w = v*u,
    x = (1-v)*u,
    y = 0.23245,
    z2 = 13,
    v1 = 0.0,
    v2 = 0.0)

function RHS(ydot,y,par,t)

    a = par.a;
    d = par[2];
    T = par.T;
    d = par[4];
    c = par.c;
    e = par.e;
    f = par.f;
    g = par.g;
    h = par.h;
    k = par.k;
    l = par.l;
    m = par.m;
    n = par.n;
    p = par.p;
    q = par.q;
    z = par[16];
    s = par.s;
    t = par.t;
    u = par.u;
    w = par.w;
    x = par.x;
    #y = par[22];
    z = par[23];
    v1= par.v1;
    v2= par.v2;

    fM = d * y[1] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
    h = c * d * y[2] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;

    ydot[1] = -fM - (v1*y[1]-v2*y[2]);
    ydot[2] = -h + (v1*y[1]-v2*y[2]);
    ydot[3] = fM - k*y[3]-(v1*y[3]-v2*y[4]);
    ydot[4] = h - k*y[4]+(v1*y[3]-v2*y[4]);
    ydot[5] = k*(y[3]-y[5])-z*y[5]-(v1*y[5]-v2*y[6]);
    ydot[6] = k*(y[4]-y[6])-z*y[6]+(v1*y[5]-v2*y[6]);
    ydot[7] = z*(y[5]+y[6])-k*y[7];
    ydot[8] = k*(y[5]-y[8])-z*y[8]-(v1*y[8]-v2*y[9]);
    ydot[9] = k*(y[6]-y[9])-z*y[9]+(v1*y[8]-v2*y[9]);
    ydot[10] = z*(y[8]+y[9])-k*y[10] + k*y[7];
    ydot[11] = k*(y[8]-y[11])-z*y[11]-(v1*y[11]-v2*y[12]);
    ydot[12] = k*(y[9]-y[12])-z*y[12]+(v1*y[11]-v2*y[12]);
    ydot[13] = z*(y[11]+y[12])-k*y[13] + k*y[10];
    ydot[14] = k*(y[11]-y[14])-z*y[14]-(v1*y[14]-v2*y[15]);
    ydot[15] = k*(y[12]-y[15])-z*y[15]+(v1*y[14]-v2*y[15]);
    ydot[16] = z*(y[14]+y[15])-k*y[16] + k*y[13];
    ydot[17] = m*k*y[14]-z*y[17]-(v1*y[17]-v2*y[18])-t*y[17];
    ydot[18] = m*k*y[15]-z*y[18]+(v1*y[17]-v2*y[18])-t*y[18];
    ydot[19] = (1-m)*k*y[14]-(z+s)*y[19]-(v1*y[19]-v2*y[20])-(p+q)*y[19];
    ydot[20] = (1-m)*k*y[15]-(z+s)*y[20]+(v1*y[19]-v2*y[20])-(p+q)*y[20];
    ydot[21] = (1-m)*k*(y[14]+y[15]);

    return ydot
end

function condition(out,y,t,integrator) # Event when event_f(u,t) == 0
    out[1] = integrator.t - integrator.p[1] - integrator.p[2];
    out[2] = integrator.p[1] - integrator.t;
end

function affect!(integrator, idx)
    if idx == 1
        integrator.p.v1 = integrator.p[8]*integrator.p[9];
        integrator.p.v2 = integrator.p[8]*(1-integrator.p[9]);
    else
        integrator.p .= 0
    end
end

cb = VectorContinuousCallback(condition,affect!,2);
y0 = [2323.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1];
tspan = (0.0, 100.0);
prob = ODEProblem(RHS,y0,tspan,par);
sol = solve(prob, Rodas5(), callback = cb, saveat = 1.0)

```

However, note that there are funny things in your model. You define `d` and `z` twice, so in the LabelledArray it caught that and threw an error. To get around that, I just renamed them `d1` and `d2`, and `z1` and `z2`, and left the way you were dereferencing in the ODE the same so you can see the oddities there (you’re grabbing it from an index for `d1`, so the other definition isn’t used? and why are some parameters like `y` non-existent?). So you probably need to clean up the definition of the ODE before it’s anything sensible, but in terms of the part for the question, this shows a clean way how to handle the parameter changes since then it’s just `integrator.p.v1 = ...` or to zero all parameters it’s just `integrator.p .= 0`.

---

<div class="post-metadata">

**Author:** ![dynova](https://avatars.discourse-cdn.com/v4/letter/d/e47c2d/32.png) [@dynova](https://discourse.julialang.org/u/dynova)\
**Post date:** [June 17, 2020, 8:21pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/3 "2020-06-17T20:21:47Z")

</div>

Thank you - that was very helpful! After I cleaned up my code, it runs properly now.

Since this was a minimal working example, I want to run my (highly) stiff system of ODEs across a wide range of parameter values. I say “highly” because of the discontinuities accounted for by the Callback functions.

As I’m running the solver in a loop over the parameter values, I noticed that the code freezes at a certain iteration (say, iteration 493 out of 100,000). This is after I supplied the full Jacobian and used Rodas5(). When I use Rodas4(), the code freezes at iteration 350.

So, my question is: I’m not sure which solver I should use, and what the values of reltol and abstol should be. Is there a general procedure to go about this?

Once again, thank you for your help!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 17, 2020, 8:24pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/4 "2020-06-17T20:24:47Z")

</div>

> [@dynova](#):
>
> As I’m running the solver in a loop over the parameter values, I noticed that the code freezes at a certain iteration (say, iteration 493 out of 100,000). This is after I supplied the full Jacobian and used Rodas5(). When I use Rodas4(), the code freezes at iteration 350.
> 
> So, my question is: I’m not sure which solver I should use, and what the values of reltol and abstol should be. Is there a general procedure to go about this?

Guess and check is usually pretty good. What do you mean by freezes? If you put an `@show t` in your function, is it getting seemingly stuck at a point? Is there something special at that point, like a discontinuity? Decreasing tolerances, like `abstol=1e-8,reltol=1e-8` can increase the stability of solving some equations so it’s worth a try.

---

<div class="post-metadata">

**Author:** ![dynova](https://avatars.discourse-cdn.com/v4/letter/d/e47c2d/32.png) [@dynova](https://discourse.julialang.org/u/dynova)\
**Post date:** [June 17, 2020, 8:31pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/5 "2020-06-17T20:31:20Z")

</div>

I see - yes, the solver is getting stuck at a point, so I’ll look into it.

---

<div class="post-metadata">

**Author:** ![dynova](https://avatars.discourse-cdn.com/v4/letter/d/e47c2d/32.png) [@dynova](https://discourse.julialang.org/u/dynova)\
**Post date:** [June 17, 2020, 9:48pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/6 "2020-06-17T21:48:27Z")

</div>

In debugging, I noticed the following:

1. “ydot” has some negative values at various time steps.
2. The ODE integrator is stepping through the full time interval that was specified, for “proper” sets of parameter values.
3. I might have specified the callbacks incorrectly…
4. I would like to view the determinant and spectral radius of the Jacobian at each iteration. How can I do this?

For #3, my code is below:

```julia
function condition(out,y,t,integrator)
    out[1] = integrator.t - integrator.p[1] - integrator.p[2];
    out[2] = integrator.p[1] - integrator.t;
end

function affect!(integrator, idx)
    if idx == 1
        integrator.p[24] = integrator.p[8]*integrator.p[9];
        integrator.p[25] = integrator.p[8]*(1.0-integrator.p[9]);
    else
        integrator.y .= 0
        integrator.jacobian .= 0
    end
end

```

I believe the issue is in the `else` clause…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 18, 2020, 2:31pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/7 "2020-06-18T14:31:40Z")

</div>

> [@dynova](#):
>
> - “ydot” has some negative values at various time steps.
> - The ODE integrator is stepping through the full time interval that was specified, for “proper” sets of parameter values.
> - I might have specified the callbacks incorrectly…

There are a few things in your equations that probably need fixing that I pointed out. Same variables are overwritten to multiple values, some aren’t used, etc. and I assume those are all accidents.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 18, 2020, 2:32pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/8 "2020-06-18T14:32:19Z")

</div>

> [@dynova](#):
>
> `integrator.jacobian .= 0`

That’s not a thing. I don’t know what documentation you’re reading to come up with that, and I don’t know what you’re trying to do there.

---

<div class="post-metadata">

**Author:** ![dynova](https://avatars.discourse-cdn.com/v4/letter/d/e47c2d/32.png) [@dynova](https://discourse.julialang.org/u/dynova)\
**Post date:** [June 18, 2020, 3:46pm UTC](https://discourse.julialang.org/t/callbacks-for-parameters-in-a-stiff-system-of-odes/41580/9 "2020-06-18T15:46:18Z")

</div>

Got it. Thanks for the feedback.

I’ve since figured out that the issue in my overall code is actually… not with the numerical integration routine. It’s elsewhere, and more manageable. It should be working soon.