DifferentialEquations with nested functions + Callbacks

New to Julia
1

Hello, now that I (sort of) know how to calculate differential equations with nested functions I would like to learn how to use it with callbacks. Callbacks are used to run a differential equation to a certain point, changing the value of a parameter and then continue running. And that’s really handy when working with seasonal data.

I know how a callback looks like when using a differential equation that does not have nested functions, but my experiments with a callback and differential equation with nested functions have failed so far.

The outline of the code (I will post the whole code if you want, but I don’t know how to make a MWE because there are many parts)

function FF (parameters)

function CalcG (parameters)

function dNdT (parameters) #change in nutrient 'biomass'
       something* CalcG (parameters)
end

function dPdT (parameters) #change in plant biomass
       something/ FF (parameters)+ CalcG (parameters)
end

function dAdT (parameters) #change in animal biomass
       function FF (parameters)
end

And then dNdT, dPdT, dAdT are placed together into one array and calculated at the same time:

function dBdT_test(dB, bioS, p, t) #where dB is change in biomass

    #extinction
    @. bioS[bioS<0.000001] = 0.0 #here I don't know if this should be placed before or after the equations
    @. dB[bioS<0.000001] = 0.0

    dN = dB[1:num_nutrient]
    dP = dB[num_nutrient+1:(arr_size-num_animal)]
    dA = dB[(arr_size-num_animal+1):arr_size]

    dB[1:num_nutrient] = dNdT(dN, D, S, N, K, v1, v2, r, bioS)
    dB[num_nutrient+1:(arr_size-num_animal)] = dPdT(dP, K, N, r, bioS, X, num_plant, num_animal, num_nutrient)
    dB[(arr_size-num_animal+1):arr_size] = dAdT(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q)

    #extinction
    @. bioS[bioS<0.000001] = 0.0 
    @. dB[bioS<0.000001] = 0.0 #here I don't know if this should be placed before or after the equations

end

tspan = (0.0,500.0)
p = (foodWeb, num_plant, num_animal, num_nutrient, D, S, N, K, v1, v2, r, X, e, bm, b, h, ω, q)

prob = ODEProblem(dBdT,bioS,tspan,p)
sol=solve(prob)

Is it good to have dNdT, dPdT, dAdT as separate functions and then call them into dBdT ?

I only want to put a callback for a specific species (lets say species 42) in dAdT.

Should I place the callback into function dAdT or function dBdT?

I tried putting it into dBdT like so, but then f1 is not defined

dB[(arr_size-num_animal+1):arr_size] = dAdT_harvest(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q) - (catch_co*bioS[42]*effort + f1)`

And the callback is

function condition(u,t,integrator)
  t == 200.
end
function affect!(integrator)
  for c in full_cache(integrator)
    integrator.p[21] = 0.01
  end
end
cb = DiscreteCallback(condition, affect!)

tspan = (0.0,500.0)

p = (foodWeb, num_plant, num_animal, num_nutrient, D, S, N, K, v1, v2, r, X, e, bm, b, h, ω, q, catch_co, effort, f1)

prob = ODEProblem(dBdT,bioS,tspan,p)
sol=solve(prob,callback=cb,tstops=[200.0])

Suggestions are greatly appreciated

2

There is no difference between putting code in functions or not, other than for aesthetics and composibility. That said, I’m not sure what your question is.

1 Like
3

Apologies for not being clear.

When I try to introduce a callback into my (working) code, Julia gives an error that f1 is not defined. I would like to know how to correctly add a callback as I can’t find a relatable example in the callback tutorial.

The code is

function dBdT_harvest(dB, bioS, p, t)

    #extinction
    @. bioS[bioS<0.000001] = 0.0
    @. dB[bioS<0.000001] = 0.0

    dN = dB[1:num_nutrient]
    dP = dB[num_nutrient+1:(arr_size-num_animal)]
    dA = dB[(arr_size-num_animal+1):arr_size]

    dB[1:num_nutrient] = dNdT(dN, D, S, N, K, v1, v2, r, bioS)
    dB[num_nutrient+1:(arr_size-num_animal)] = dPdT(dP, K, N, r, bioS, X, num_plant, num_animal, num_nutrient)
    dB[(arr_size-num_animal+1):arr_size] = (dAdT(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q) - (catch_co*bioS[42]*effort + f1)) # I introduce f1 here 

    #extinction
    @. bioS[bioS<0.000001] = 0.0
    @. dB[bioS<0.000001] = 0.0

end

function condition(u,t,integrator)
  t == 200.
end

function affect!(integrator)
  for c in full_cache(integrator) # I do not understand the purpose of this line
    integrator.p[21] = 0.01 # Here f1 is the 21st parameter
  end
end

cb = DiscreteCallback(condition, affect!)

tspan = (0.0,500.0)

p = (foodWeb, num_plant, num_animal, num_nutrient, D, S, N, K, v1, v2, r, X, e, bm, b, h, ω, q, catch_co, effort, f1) # the error happens here

prob = ODEProblem(dBdT_harvest,bioS,tspan,p)
sol=solve(prob,callback=cb,tstops=[200.0])

Error is UndefVarError: f1 not defined in top-level scope at base\none. I assume the issue is something with this part:

function affect!(integrator)
  for c in full_cache(integrator) 
    integrator.p[21] = 0.01
  end
end
4 
function affect!(integrator)
  integrator.p[21] = 0.01
end

but that’s not the line causing your issue.

dB[(arr_size-num_animal+1):arr_size] = (dAdT(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q) - (catch_co*bioS[42]*effort + f1)) # I introduce f1 here

but f1 isn’t defined, so :man_shrugging:

5

I thought f1 corresponds to

function condition(u,t,integrator)
  t == 200.
end

function affect!(integrator)
     integrator.p[21] = 0.01 # Here f1 is the 21st parameter
end

My goal is to change the value of catch_co to 0.01 at the 200th time step.

6

You never said f1 = ... I don’t even know what f1 is. All of dN, D, S, N, K, v1, v2, r are also not defined in your code. Are they from p? You need to write down any namings you come up with, otherwise both me and the computer will have no idea f1 = p[21].

7

I am trying to follow the examples here and here so I thought it was common to use f1.

Everything in the code is defined and works, except for f1 which I now introduce as part of the callback. I will think about this some more and try to make my question more understandable.

8

f1 is a field in the array type in that example. You’re using normal arrays so it doesn’t exist.

9

Okay, I got it to run.

For people struggling with this,

First, make sure the parameter (not the thing/variable the differential equation is supposed to calculate) you are going to change is defined somewhere already, even if the starting value is 0)

For example, catch_co = 0

Then for the call back

function condition(u,t,integrator)
  t == 200.
end

Where 200 is the time step you want your parameter to change

Then you have to specify what parameter is changing and what it is changing to

function affect!(integrator)
  for c in full_cache(integrator)
    integrator.p[19] = 0.01
  end
end

Where [19] corresponds to the number of your parameter in your p list (see below), and it changes to a value of 0.01.

cb = DiscreteCallback(condition, affect!)

tspan = (0.0,500.0)

p = (foodWeb, num_plant, num_animal, num_nutrient, D, S, N, K, v1, v2, r, X, e, bm, b, h, ω, q, catch_co, effort) 
#In this case `catch_co = 0` is the 19th parameter in the list.

prob = ODEProblem(dBdT_harvest,bioS,tspan,p)
sol=solve(prob,callback=cb,tstops=[200.0])

So that’s how a callback with nested functions is done.