Returning a value calculated inside the ODE function

Specific Domains Modelling & Simulations
1

Hi

I am using julia via R to simulate compartmental epidemic ODE models. Is there a way to return a value that is calculated inside the ODE function for every step of the integration? So far it only returns the states. Below is a toy model, and I would like to return the value of “FOI” and “N” at each time step in addition to the states. What is the best way to do this? return (du[1], du[2], du[3], FOI, N) does not work

library(JuliaCall)
library(diffeqr)

func_julia2 <- JuliaCall::julia_eval("
function func_julia2(du, u , p , t)
  
  beta = p[1]
  gamma = p[2]
  mu = p[3]
  sigma = p[4]
  eta = p[5]
  phi = p[6]
  
  S = u[1]
  I = u[2]  
  R = u[3]
  N = S + I + R
  
  beta_eff = beta * (1.0+eta*sin(pi*(t-phi)/365.0)^2.0)
  FOI = beta_eff*I/N
  
  du[1] = -FOI*S + mu*N - mu*S + sigma*R
  du[2] = FOI*S - gamma*I - mu*I
  du[3] = gamma*I - mu*R - sigma*R

end")

JuliaCall::julia_assign("u", c(10000,1,0))
JuliaCall::julia_assign("p", c(0.25, 1/5, 1/(70*365), 1/365, 0.2, 10))
JuliaCall::julia_assign("t", c(0.0, 3650.0))
prob3 = JuliaCall::julia_eval("ODEProblem(func_julia2, u, t, p)")
sol = enviro$solve(prob3, enviro$AutoVern7(enviro$KenCarp4()), 
                   saveat=1.0, abstol=1e-8, reltol=1e-6)
2

There is a SavingCallback documented here Callback Library · DifferentialEquations.jl which can do what you want. I find it a little awkward to use in that the saved variables aren’t in the same place as the state variables and you don’t get interpolation. Sometimes it’s easier to make helper functions that recreate things from the state variables, and in cases where I’m using a callback to simulate a control system I’ll add dummy variables to the states (with the rate of change set to zero) and update their values from inside the callback.

1 Like