Is there a way to run differential equations to a specific point and then change input values, continue running, and graph it together?

New to Julia
,
1

I have a simple 3-species food chain that I would like to run for 200 steps so the system stabilizes. I would then like to start harvesting and show how populations fluctuate. Is there a way to do this?

The natural/non-harvested system is:

using DifferentialEquations
f = @ode_def_bare Pristine begin
dBB = r*(1- BB/K)BB -xIyIBBI((BB^h/(B0^h + BB^h))/eIB)
dBI = -xIBI + xIyIBBI(BB^h/(B0^h + BB^h))-xT yTIBT*((BI^h/(B0^h + BI^h))/eTI)
dBT = -xTBT + xTyTIBT(BI^h/(B0^h + BI^h))-qBTE
end K r yIB yTI eIB eTI B0 h xI xT q E
u0 = [155.0,107.0,93.0]
tspan = (0.0,1000.0)
p = (700, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0, 0)
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23())

And it looks like this (sorry at the moment I haven’t figured out how to add a main title to the plots and not for each subplot):

image
image841×564 33.6 KB

Plot code is

using Plots
plot(sol,xlabel = “Time” ,ylabel = “Density”, title = “Food Chain Pristine”, lw=0.5, layout = (3,1))

Then I would like to start harvesting after 200 steps so q = 0.01 and E = 50 (instead of both 0 in the pristine scenario).

Is this possible?

2

The easiest way is to just make it a single diffeq solve with a callback that implements the change.

3

Thank you for your response here and on my user_cache question. I have tried to implement the callback using the Callback example but I don’t know how to have two parameters changed. Here is my attempt which generates ERROR: LoadError: invalid redefinition of constant f.

I probably need to put an f2 component here

using DifferentialEquations
mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T

and another error source is probably stating initial values of q and E as 0.0 while specifying 0.01 and 50 under parameters. q and E do not change over time in this model (they will in future ones) so perhaps they should not be put into u0?

u0 = [155.0,107.0,93.0,0.0,0.0]
tspan = (0.0,1000.0)
p = (450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50)

Whole code:

using DifferentialEquations
mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end
f = @ode_def_bare ManagedLinear begin
dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E + u.f1 +u.f2
end K r yIB yTI eIB eTI B0 h xI xT q E
const tstop1 = [200.]

function condition(u,t,integrator)
  t in tstop1
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 0.1
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f2 = 50
  end
end
save_positions = (true,true)

cb = DiscreteCallback(condition, affect!, save_positions=save_positions)

save_positions = (false,true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)
u0 = [155.0,107.0,93.0,0.0,0.0]
tspan = (0.0,1000.0)
p = (450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50)
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23())

using Plots

plot(sol, title = "Food Chain Managed Linear", xlabel = "Time" ,ylabel = "Density", lw=0.5, layout = (3,1))
4

That happens when you define a function twice, the first time having it be a constant. Try your code in a fresh REPL session.

For your code, your issues are you haven’t defined condition2 (which I assume needs a second tstop and f2 handling), and you need to add the callback to the solving: sol=solve(prob,Rosenbrock23(),callback=cbs).

Note that you can do this without DEDataVectors though. This is probably an example of how you want to work:

using DifferentialEquations

f = @ode_def_bare ManagedLinear begin
  dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
  dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
  dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E + f1
end K r yIB yTI eIB eTI B0 h xI xT q E f1

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

u0 = [155.0,107.0,93.0]
tspan = (0.0,1000.0)
p = [450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50, 0.0] # initial f1 = 0.0
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23(),callback=cb)

using Plots

plot(sol, title = "Food Chain Managed Linear", xlabel = "Time" ,ylabel = "Density", lw=0.5)
5

I have included a condition2, ‘affect2’ and grouped it in a ‘CallBackSet’, but when I run it, the plot generated does not differ from the original code without the callback. What am I missing?

using DifferentialEquations

f = @ode_def_bare ManagedLinear begin
  dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
  dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
  dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E + f1 + f2
end K r yIB yTI eIB eTI B0 h xI xT q E f1 f2

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

function condition2(u,t,integrator)
  t == 200.
end
function affect2!(integrator)
  for c in full_cache(integrator)
    integrator.p[14] = 50
  end
end
cb2 = DiscreteCallback(condition2, affect2!)

cbs = CallbackSet(cb,cb2)

u0 = [500.0,200.0,100.0]
tspan = (0.0,1000.0)
p = [450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50, 0, 0] # initial f1 = 0.0 f2 = 0.0
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23(),callbacks=cbs)

using Plots

plot(sol, title = "Food Chain Managed Linear Callback", xlabel = "Time" ,ylabel = "Density", lw=0.5)
6

My bad here. I forgot to set the tstops so that it hits t=200.0 exactly to trigger the callbacks.

sol=solve(prob,Rosenbrock23(),callbacks=cbs,tstops=[200.0])

Hope that makes sense.

7

Hello again,

Even with the tstops it generates the same graph. Perhaps the problem is with…

p = [450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50, 0, 0] # initial f1 = 0.0 f2 = 0.0
# K r yIB yTI eIB eTI B0 h xI xT q E f1 f2

…since q and E should both be 0 until t=200.0 and then become 0.01 and 50?

When I set them both to 0 it generates the Pristine graph, which is what I want, but then nothing happens once t=200.0.

using DifferentialEquations

f = @ode_def_bare ManagedLinear begin
  dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
  dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
  dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E + f1 + f2
end K r yIB yTI eIB eTI B0 h xI xT q E f1 f2

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

function condition2(u,t,integrator)
  t == 200.
end
function affect2!(integrator)
  for c in full_cache(integrator)
    integrator.p[14] = 50
  end
end
cb2 = DiscreteCallback(condition2, affect2!)

cbs = CallbackSet(cb,cb2)

u0 = [500.0,200.0,100.0]
tspan = (0.0,1000.0)
p = [450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50, 0, 0] # initial f1 = 0.0 f2 = 0.0
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23(),callbacks=cbs,tstops=[200.0])

using Plots

plot(sol, title = "Food Chain Managed Linear ", xlabel = "Time" ,ylabel = "Density", lw=0.5)
8

The issue is that you used the argument callbacks instead of callback. That was my bad (and we really need to error for that). A corrected code is:

using DifferentialEquations

f = @ode_def_bare ManagedLinear begin
  dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
  dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
  dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E + f1 + f2
end K r yIB yTI eIB eTI B0 h xI xT q E f1 f2

function condition(u,t,integrator)
  t == 200.
end
function affect!(integrator)
  integrator.p[13] = 0.01
end
cb = DiscreteCallback(condition, affect!)

function condition2(u,t,integrator)
  t == 200.
end
function affect2!(integrator)
  integrator.p[14] = 50
end
cb2 = DiscreteCallback(condition2, affect2!)

cbs = CallbackSet(cb,cb2)

u0 = [500.0,200.0,100.0]
tspan = (0.0,1000.0)
p = [450, 1.1, 10, 10, 0.66, 0.85, 80, 1.2, 0.15, 0.06, 0.01, 50, 0, 0] # initial f1 = 0.0 f2 = 0.0
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,Rosenbrock23(),callback=cbs,tstops=[200.0],reltol=1e-6)

using Plots

plot(sol, title = "Food Chain Managed Linear ", xlabel = "Time" ,ylabel = "Density", lw=0.5)

plot