# Track variables within DifferentialEquations.jl function

**URL:** <https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777>\
**Category:** New to Julia\
**Created:** [March 9, 2020, 10:57pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777 "2020-03-09T22:57:40Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 9, 2020, 10:57pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/1 "2020-03-09T22:57:40Z")

</div>

Hi,

I have a simple function to solve:

```julia
using DifferejntialEquations, Plots, BenchmarkTools

var = []
function coupledd!(du, u, p, t)
    α, β, D1, D2 = p
    c1, c2 = u
    dc1 = α*c1 + D1 * (c2-c1)
    dc2 = β*c2 + D2 * (c1-c2)
    push!(var, c2)
    du .= [dc1, dc2] # note the broadcast operator!
end

u0 = rand(2)
tspan = (0.,10.)
pars = [-.2, -.1, 1., 1.]
prob = ODEProblem(coupledd!, u0, tspan, pars)
sol = solve(prob, Tsit5())
println("size(var)= ", size(var))
println("nb time steps= ", size(sol.t))

```

size(var) returns 696  
size(sol.t) return 12

However, when I print var, I get the output:  
87-element Array{Any,1}

which seems to imply that Tsit5() calls the right-hand side 8 times every time step.

What I am interested in is the value of a variable at every time step (once per time step). So I would like to monitor variables. Is there a facility for that? Thanks.

Gordon

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 10, 2020, 1:23am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/2 "2020-03-10T01:23:15Z")

</div>

You can either specify the keyword `saveat`, or just call the solution using the interpolator, e.g.

```julia
sol(t) # Gives you the solution at time `t` through interpolation of the same order as the solution method

```

instead of `sol.t`

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 1:28am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/3 "2020-03-10T01:28:25Z")

</div>

I was not clear. My mistake for choosing to monitor the actual solution.

Say I defined the variable z=c1+c2 in the coupledd!!() function and wish to store z in a vector at the same time levels as sol.t . Of course, I could recompute z outside the function, but why should I, especially if z is expensive to compute. I just asked a minimal example.

Thanks.

---

<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:** [March 10, 2020, 5:47am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/4 "2020-03-10T05:47:47Z")

</div>

> [@erlebach](#):
>
> Say I defined the variable z=c1+c2 in the coupledd!!() function and wish to store z in a vector at the same time levels as sol.t . Of course, I could recompute z outside the function, but why should I, especially if z is expensive to compute. I just asked a minimal example.

You can use a `DiscreteCallback` for this, and have it fire after every step by making it `true`. There’s actually already a premade callback in the library for this exact case, the `SavingCallback`.

[https://docs.juliadiffeq.org/latest/features/callback\_library/#saving\_callback-1](https://docs.juliadiffeq.org/latest/features/callback_library/#saving_callback-1)

Or you can use the integrator interface for stepwise control of the solution:

[https://docs.juliadiffeq.org/latest/basics/integrator/](https://docs.juliadiffeq.org/latest/basics/integrator/)

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 10:00am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/5 "2020-03-10T10:00:39Z")

</div>

Thank you, Chris. My function is much more complicated. But it does appear that I must recompute the variables of interest, and redefine them in a new function in order to do the monitoring. That violates the DRY principle, and can be error-prone. It would be nice to tag the variables of interest in the function defining the right hand side and have the callbacks generated automatically. Am I making sense?

Thanks.

---

<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:** [March 10, 2020, 11:15am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/6 "2020-03-10T11:15:06Z")

</div>

> [@erlebach](#):
>
> But it does appear that I must recompute the variables of interest, and redefine them in a new function in order to do the monitoring

Not necessarily. You’ll only be re-evaluating at the final `u` if the method is FSAL, which many methods are not (though B-stable methods have it). So the callback, i.e. the value at t + dt, will be a brand new `u`. You could then Ref these over to `p` if you really wanted to. I can throw an example together if you think it’s worth it.

But I am not sure it is. You would decrease your calculations by at most 1/8 if this is Tsit5 for example. If this is something for stuff equations, this would be even less. For example, with Rodas5 this would be saving 1/(N+9) where N is the number of states in the system (or if sparse, `maximum(colorvec)`). I think the case where this factor is largest is BS3 where it is an FSAL algorithm but 4 stages, so you’d reduce the calculation of the explicit algebraic variables by 1/5. So in a worst case scenario, where the calculation of the explicit algebraic variables is 100% of your code’s time and you use BS3, then you at most have a 20% slowdown from this. In total, there are probably a million other things to optimize first.

Is there an interface that you’re used to that allows/tracks explicit algebraic variables in a nice way?

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 12:14pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/7 "2020-03-10T12:14:26Z")

</div>

Brian2 is such an interface, using Python. I may have to provide an example. Your idea of using p is what I need I think, if by p you mean the parameter array. I am using these variables for debugging purposes with visualization. It is not a out thecomputer time, but about not wishing to introduce new errors.i do not wish to defi e complex e pressions undergoing structural change in two different locations.

Thanks.

Gordon

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 4:10pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/8 "2020-03-10T16:10:23Z")

</div>

Chris, here are some more details.

```julia
function myCallback(u, t, integrator)
	C = @view(u[1:1:2]) # Calcium
    Ce = @view(u[3:1:4]) # Calcium in ER
    I = @view(u[5:1:6]) # IP3
    h = @view(u[7:1:8])
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
end

```

I am returning several variables. There are 8 equations, two equations for C, two for Ce, etc.

I setup the callbacks:

```julia
saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)

prob = ODEProblem(Cresswell!, u0, tspan, pars) # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)

```

I thought that Vector::Any was completely general in the sense that Vector elements could be any type. But that is not the case according to the first line of the error message. Here is a zoomed view:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/a/bafcd413219df45ebc33972a3cf17e32a4fff4d0.png)

Each of the return elements is an array of two elements (could be more in general). I want the capability of adding more variables without having to change the return type. So I am concerned with usability, not efficiency, since this is a debugging tool for me at this stage.

I did check out [https://docs.juliadiffeq.org/release-5.0/features/callback\_library.html](https://docs.juliadiffeq.org/release-5.0/features/callback_library.html)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/4/444bcf1a3c0db3bc0ca829e2610431d98f6d9589.png)  
(I cannot copy/paste this error. The part in read is not copiable. Why would that be?)

I would like to also know how bring in my parameter array “p” into the myCallback() function. I could access it as a global variable, but that is poor programming.

Thanks!

* * *

Here is the full code.

```julia
using DifferentialEquations, Plots, BenchmarkTools

function globalParametersEvan()
    # Try to duplicate Evan's model as much as possible
    d = Dict()
    d[:Rgas] = 8.31 # joule / kelvin / mole
    d[:D_I] = 5.3e-10 # Hz * meter^2
    d[:D_C] = 5.3e-10 # Hz * meter^2
    d[:D_CE] = 5.3e-10 # Hz * meter^2
    d[:Cm] = 1.e-6 # Farad * meter^-2
    d[:e] = 1.602176634e-19 # Amp * second
    d[:N_A] = 6.02214076e23 # mole^-1
    d[:k_B] = 1.38064852e-23 # meter^2 * Kg * Hz^2 / kelvin
    d[:F] = 96485.3329 # Amp * second / mole = Cb / mole
    d[:V_T] = 0.026 # voltage_cubic
    d[:P_Ca] = 4.46e-15 # meter * Hz
    d[:Crest] = 0.073e-6 # Molar
    d[:Vrest] = -80. # Volt
    d[:Ce] = 1.800 # mole/meter**3]
    d[:Ce] = 1.e-4 # mole/meter**3] # WHICH ONE TO CHOOSE?

    d[:O_delta] = 0.6e-3 # Hz * mole/meter^3
    d[:k_delta] = 1.5e-3 # mole/meter^3
    d[:K_delta] = 0.1e-3 # mole/meter^3
    d[:K_D] = 0.7e-3 # mole/meter^3
    d[:K_3] = 1.0e-3 # mole/meter^3
    d[:O_3K] = 4.5e-3 # Hz*mole/meter^3
    d[:o_5P] = 0.05 # Hz
    d[:K_5P] = 10.0e-3 # mole/meter^3
    d[:o_delta] = 0.6e-3 # mole/meter^3
    d[:k_delta] = 1.5e-3 # Hz * mole/meter^3
    d[:K_delta] = 0.1e-3 # mole/meter^3
    d[:Lambda] = 2100e-18 # meter^3
    d[:rho_A] = 0.18 # ER-to-cytoplasm volume ratio ?
    d[:rho] = d[:rho_A] / (1. + d[:rho_A]) # # ER-to-cytoplasm volume ratio ?

    # Multiply value of amplitudes by cytosolic volume (changed rho_A to rho)
    # It does not look as if I did the multiplication
    #self.o_delta = 0.6 # * umolar * Lambda * (1-rho) / second
    d[:o_delta] = 0.6e-3 # * mole/meter**3 * Lambda * (1-rho) / second
    d[:o_3K] = 4.5e-3 # * mole/meter**3 * Lambda * (1-rho) / second

    d[:Omega_L] = 0.1 #/second # Maximal rate of Ca^2+ leak from the ER
    d[:Omega_2] = 0.2 #/second # IP_3R binding rate for Ca^2+ inhibition
    d[:Omega_u] = 0.2-3 #/second # uptake/dissassociation constant

    # --- IP_3R kinectics
    d[:d_1] = 0.13e-3 #*mole/meter**3 # IP_3 binding affinity
    d[:d_2] = 1.05e-3 #*mole/meter**3 # Ca^2+ inactivation dissociation constant
    d[:O_2] = 0.2e3 #/(mole/meter**3)/second # IP_3R binding rate for Ca^2+ inhibition
    d[:d_3] = 0.9434e-3 #*mole/meter**3 # IP_3 dissociation constant
    d[:d_5] = 0.08 #*mole/meter**3 # Ca^2+ activation dissociation constant

    d[:p_open] = 1
    d[:P_r] = 1e-3 #*mole/meter**3/second
    d[:P_CE] = 1e-3 #*mole/meter**3/second

    d[:eta_p] = 1
    d[:d_ER] = 1e-6 #* meter
    d[:dv_ER] = 10e-3
    d[:K_P] = 0.05e-3

    d[:L] = 8e-6 # umeter
    d[:R] = 3e-6 # meter
    d[:r] = 2e-6 # meter
    d[:O_delta] = 0.6e-21 # mole/meter**3 * Lambda * (1-rho) / second
    d[:O_5P] = 0.05e-21 # mole/meter**3 * Lambda * (1-rho) / second
    d[:F] = 0.1 # ALREADY DEFINED?
    d[:D] = 0.05

	# Evan parameters
	d[:vIPR] = 6 # sec^-1
	d[:vSERCA] = 4.4e-3 # mole m^-3 s^-1 (1 Molar = 1 Mole/liter = 10^3 mole/meter^3)
	d[:vLeak] = 0.11 # sec^-1
	d[:d1] = 0.13e-3 # mole m^-3
	d[:d5] = 0.08234e-3 # mole m^-3
	# open channels
	d[:d2] = 1.049e-3 # mole m^-3
	d[:d1] = 0.13e-3 # mole m^-3
	d[:d3] = 0.9434e-3 # mole m^-3
	d[:a2] = 0.20e-3 # mole m^-3 s^-1
	# IP3 (J_δ)
	d[:o_δ] = .025e-3 # Mole m^-3 s^-1
	d[:K_δ] = 0.5e-3 # Mole m^-3, Ca2+ affinity of PLCδ
	d[:k_δ] = 1.e-3 # Mole m^-3, IP3 inhibiting affinity // Oschmann2009: 1.5 | COMPGLIO pg130: 1.
	# End Evan Parameters
    return d
end

lJp = []
lJr = []
lJ1 = []
lJleak = []
lJSERCA = []
lJIPR = []
lJβ = []
lCount = []
lTime = []

function getCurrents(C, Ce, I, h, t, p)
	minf = @. Hill(I,p[:d1],1) * Hill(C,p[:d5],1)
	Jleak = @. p[:vLeak] * (Ce-C)
	JIPR = @. p[:vIPR] * (minf*h)^3 * (C-Ce)
	JSERCA = @. p[:vSERCA] * Hill(C,p[:d5],2)
	J_β = sin(2*π*t)^10 #*mole/meter **3/second : mole/meter** 3/second
	J_δ = @. p[:o_δ] * Hill(p[:k_δ], I, 1) * Hill(C, p[:K_δ], 2)
	J5P = @. p[:o_5P] * Hill(I, p[:K_5P], 1) #: mole/second # o_5P, K_5P
	J3K = @. p[:o_3K] * Hill(C, p[:K_D], 4) * Hill(I, p[:K_3], 1) #: mole/second # o_3K, K_D, K_3
	Q2 = @. p[:d2] * (I + p[:d1]) / (I + p[:d3])
	OmegaH = @. p[:a2] * (Q2 + C)
	hinf = @. p[:d_2] * (I + p[:d_1]) / (p[:d_2]*(I + p[:d_1]) + (I + p[:d_3])*C) #: 1 # d_2, d_1, d_3
	return J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2
end

function Cresswell!(du, u, p, t)
    C = @view(u[1:1:2]) # Calcium
    Ce = @view(u[3:1:4]) # Calcium in ER
    I = @view(u[5:1:6]) # IP3
    h = @view(u[7:1:8]) # fraction of open channels
    push!(lCount, 1)

    # rhs_C, etc. Are arrays, created via implicit vectorization involving numpy arrays
    dR2 = p[:R]^2 - p[:r]^2 #: meter**2
    s = p[:R] + p[:r]
    V = @. p[:Vrest] + (p[:F]*s/p[:Cm]) * (C - p[:Crest]) #: volt
    VVT = -2. * V/p[:V_T] #: 1 # value of 0, which leads to a zero denominator

	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, p)
 	p[:J_δ] = J_δ
    push!(lJleak, Jleak)
    push!(lJIPR, JIPR)
    push!(lJSERCA, JSERCA)
    push!(lJβ, J_β)
	push!(lTime, t)

    dh = @. OmegaH * (hinf - h) #: 1
    dI = @. (J_β + J_δ - J3K - J5P)

    # GE: divde Jp by Lambda (2020-03-09)
    dC = @. (JIPR - JSERCA + Jleak) #: mole / meter**3
    dCE = @. -(1. / p[:rho_A]) * (JIPR - JSERCA + Jleak) #: mole / meter**3
    dI = @. (1. *J_β + 1. * J_δ - 1. *J3K - 1. *J5P) # / (p[:Lambda]*(1-p[:rho])) #: mole/meter**3
    dh = @. OmegaH * (hinf - h) #: 1

    du[1:2] = dC
    du[3:4] = dCE #[0,0]
    du[5:6] = dI # [0,0] # blows up!
    du[7:8] = dh
end

nb_neurons = 2
eqs_per_neuron = 4
# Initial Conditions
function initialConditions()
    c0 = rand(nb_neurons)
    ce0 = rand(nb_neurons)
    h0 = rand(nb_neurons)
    I0 = rand(nb_neurons)
    # Initial condition vector
    u0 = vcat(c0,ce0,h0,I0)
end

function callback(u, t, integrator)
	C = @view(u[1:1:2]) # Calcium
    Ce = @view(u[3:1:4]) # Calcium in ER
    I = @view(u[5:1:6]) # IP3
    h = @view(u[7:1:8])
	# How to I get my parameter array as a local variable? 
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
end

getParams = globalParameters_mod.globalParametersEvan

saved_values = SavedValues(Float64, Vector::Any)
cb = SavingCallback(callback, saved_values)

u0 = initialConditions()
pars = getParams()
tspan = (0., 300.) # make these real
#u0 = [.1, .2, .3, .4]``
# We have 8 equations. How to collect them?
prob = ODEProblem(Cresswell!, u0, tspan, pars) # d is a parameter dictionary
sol = solve(prob, Tsit5(), callback=cb)

```

Here is a modified version of the callback unction:

```julia
function callback(u, t, integrator)
	C = @view(u[1:1:2]) # Calcium
    Ce = @view(u[3:1:4]) # Calcium in ER
    I = @view(u[5:1:6]) # IP3
    h = @view(u[7:1:8])
	# How to I get my parameter array as a local variable?
	J5P, J3K, J_β, J_δ, Jleak, JIPR, JSERCA, hinf, OmegaH, minf, Q2 = getCurrents(C, Ce, I, h, t, pars)
	cur = Dict()
	cur[:J5P] = J5P
	cur[:J3K] = J3K
	cur[:J_β] = J_β
	cur[:J_δ] = J_δ
	cur[:Jleak] = Jleak
	cur[:JIPR] = JIPR
	cur[:JSERCA] = JSERCA
	cur[:hinf] = hinf
	cur[:OmegaH] = OmegaH
	cur[:minf] = minf
	cur[:Q2] = Q2
	return cur
end

saved_values = SavedValues(Float64, Dict{Any,Any})
cb = SavingCallback(callback, saved_values)

```

and this “seems” to work. I do not yet know how to use this callback once the solver has completed.  
But my immediate comment is that my callback is now substantially longer than it was and more unwieldy. I could have populated a dictionary in the getCurrents() method, but then I would have to access these dictionary values in the ODE solver. I assume that accessing dictionary values is very cheap?

I have another question. If I had an array:

```julia
vars = [:J5P, :J3K]

```

for example, would it be possible to create a macro to execute:

cur[:J5P] = J5P  
cur[:J3K] = J3K

?

Thanks,

---

<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:** [March 10, 2020, 6:37pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/9 "2020-03-10T18:37:29Z")

</div>

> [@erlebach](#):
>
> saved\_values = SavedValues(Float64, Vector::Any)

It’s `Vector{Any}`

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 7:06pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/10 "2020-03-10T19:06:57Z")

</div>

I was close. Yes, it works. I still do not grok Julia types. But so far so good. Now, how to use the results from the callback.

---

<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:** [March 10, 2020, 7:12pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/11 "2020-03-10T19:12:09Z")

</div>

`T{X}` is the family of `T` “of” `X`. I.e. Vector of Float64 is `Vector{Float64}`. `::` is used for declarations, like `x::T` i.e. “x is a T”.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 10, 2020, 10:16pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/12 "2020-03-10T22:16:43Z")

</div>

So in your simple example case, you have:

\frac{dc\_1}{dt} = \alpha c\_1 + D\_1(c\_2-c\_1) \\ \frac{dc\_2}{dt} = \beta c\_2 + D\_2(c\_1-c\_2)

And you want to compute, say z = c\_1 + c\_2, in addition, with the same time resolution as for c\_1 and c\_2? Without having to recompute z in post processing?

Why not pose the problem as a DAE? So you define an error vector e as:

e\_1 = -\frac{dc\_1}{dt} + \alpha c\_1 + D\_1(c\_2-c\_1) \\ e\_2 = - \frac{dc\_2}{dt} + \beta c\_2 + D\_2(c\_1-c\_2) \\ e\_3 = -z + c\_1 + c\_2

and aim to make the errors e\_1, e\_2, e\_3 zero.

The following code does this (I use `err` instead of `out` — `out` is used in the documentation of DifferentialEquations.jl):

```julia
using DifferentialEquations, Plots, BenchmarkTools
using Sundials # package of DAE solver(s)

#var = []
function coupledd!(err,du, u, p, t)
    α, β, D1, D2 = p
    #c1, c2 = u
    err[1] = -du[1] + α*u[1] + D1 * (u[2]-u[1])
    err[2] = -du[2] + β*u[2] + D2 * (u[1]-u[2])
    err[3] = -u[3] + u[1] + u[2]
# push!(var, c2)
    #du .= [dc1, dc2] # note the broadcast operator!
end

u0 = rand(3)
u0[3] = u0[1] + u0[2] # making sure that the initial values u0 satisfies the equations
du0 = zeros(3) # Not bothering about checking that du0 is correct or not
d_vars = [true, true, false] # Specifying which of the u variables are differentiable
tspan = (0.,10.)
pars = [-.2, -.1, 1., 1.]
prob = DAEProblem(coupledd!, du0, u0, tspan, pars,differential_vars=d_vars)
sol = solve(prob,IDA());

```

You can finally plot the result:

```julia
plot(sol,vars=[(0,1),(0,2),(0,3)])

```

which looks as follows:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/f/ffa75a8b35aa7813683961dd57b4698bf987a2cb.png)

Is this what you wanted to achieve?

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 10, 2020, 10:56pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/13 "2020-03-10T22:56:43Z")

</div>

For this simple example, you are correct. But I posted the actual code a few messages back. The DAE approach is not really what I want. My expression z is in reality auxiliary variables used to define du. I simply would like to easuly specify a list of variables to monitor, with an expression such as:

Mon = monitor([t, var1, var2, var3])

Then, solve the equation with solve(), and then I should be to plot using

Plot(Mon.t, Mon.var1)

There seems to be nothing around that can do this. To add a variable to monitor, simply add a variable to the list. Brian2 has this capability in Python.

---

<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:** [March 11, 2020, 2:15am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/14 "2020-03-11T02:15:39Z")

</div>

> [@erlebach](#):
>
> There seems to be nothing around that can do this. To add a variable to monitor, simply add a variable to the list. Brian2 has this capability in Python.

The ODE solvers that Brian2 uses (GSL) don’t have this capability either, so whatever workaround you’re doing is similar to what Brian2 is doing. Going to a DAE might be overkill unless the problem is already stiff (since then you have to use an implicit method), but putting the algebraic part inside of a function and just calling that function in two spots is DRY without performance issues (and most likely exactly what they are doing).

There is an issue to investigate adding this to the solvers, but there is are good reasons for not supporting it directly in the solvers too, and leaving that as a user-level detail.

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 11, 2020, 2:34am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/15 "2020-03-11T02:34:11Z")

</div>

I could imagine adding a line inside the rhs routine stating the names of the variables to be monitored along with an indication ad to whether they should be tracked on entry or on exit. In Brian, the equations are written inside a string that is parsed, and so there is more flexibility. The code is translated to Python or C.

The callback mechanism is good for the examples shown, but not for my use case. They are unwieldy when monitoring many variables and it is more difficult than it could be. The notation Mon.var returning a time series vector is very useful.

I enjoyed the discussion and will think some more.

Gordon

---

<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:** [March 11, 2020, 3:26am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/16 "2020-03-11T03:26:48Z")

</div>

> [@erlebach](#):
>
> I could imagine adding a line inside the rhs routine stating the names of the variables to be monitored along with an indication ad to whether they should be tracked on entry or on exit. In Brian, the equations are written inside a string that is parsed, and so there is more flexibility. The code is translated to Python or C.

That seems to me like a modeling layer on top of the differential equation solver (like how Brian2 has it), rather than something inside of a differential equation solver itself. It sounds like something that could have a million variations for every discipline to do something differently. Also anything with “tracking” requires careful compile time work or is slow, so that’s very likely not going to be in the straight solvers which in theory should be something that can compile down to Float64 code to deploy to other languages.

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 11, 2020, 10:21am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/17 "2020-03-11T10:21:15Z")

</div>

You make a good point, Chris . However, a good the beauty of ajulia should combin flexibility and ease of use. When my kind of callbacks or monitoring are used, for debugging, ease of use is Paramount. In other cases, efficiency is paramount. Thanks for the insight!

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [March 11, 2020, 10:33am UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/18 "2020-03-11T10:33:18Z")

</div>

I think I can get closer to what I want by using integrate seperately, and integrate step by step, with fixed time step, for debugging, and after every step, place the monitored variables in some type of dataframe.

I just came across the following blog that provides a high level overview of tabular data, which might help me: [DataFrames Archives | juliabloggers.com](https://www.juliabloggers.com/tag/dataframes/)

---

<div class="post-metadata">

**Author:** ![Luis\_Alejandro\_Marqu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luis_alejandro_marqu/32/35607_2.png) [@Luis\_Alejandro\_Marqu](https://discourse.julialang.org/u/Luis_Alejandro_Marqu)\
**Post date:** [April 21, 2022, 11:11pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/19 "2022-04-21T23:11:56Z")

</div>

Maybe a few years too late, but I found the same problem. My solution was this

```julia
using DifferentialEquations
tspan = (0.0,10.0)
dt = 0.1
sampletimes = tspan[1]:dt:tspan[2]

varz = []
function f(du,u,p,t)
    z = -2*u[1]
    if t in sampletimes
	push!(varz,z)
    end
    du[1] = z
end
u0 = [10.0]
prob = ODEProblem(f,u0,tspan,p=nothing;saveat=dt,tstops=sampletimes)
sol = solve(prob,Tsit5())

#this removes the data at duplicate times, even though I do not know why it appears
k = Int(size(varz)[1]/size(sampletimes)[1])
varz = varz[k:k:size(varz)[1]]

```

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 12, 2022, 11:06pm UTC](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777/20 "2022-07-12T23:06:22Z")

</div>

Hi Luis. Better late than never.  
I have got a similar problem, I’m a new Julia user but matlab expert(I am tired of paying them).

I want to call back certain values on my ode as follows from matlab

```julia
function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=20;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);
vel(1)=u;
vel(end)=0;
dudx=[(vel(2)-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-vel(end-1))/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+(0.05*(dudx)).^2).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);
eqn(N+1)=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

    
end

```

In matlab, we just have to put the variables as outputs after the eqn, then write a callback function (a bit complicated version of my own)

```julia
[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol) navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);
Force=cell2mat(Force);
Fr=cell2mat(Fr);

```

where as it gives me back a cell of the values with the time points, and then convert to a matlab simple array.

Can you help me achieve this in Julia?

[Next page](https://discourse.julialang.org/t/track-variables-within-differentialequations-jl-function/35777.md?page=2)