R/dyn_solve.R

In dynemu: Emulation of Dynamic Simulators via One-Step-Ahead Approach

#' Solving ordinary differntial equations (ODEs) on a specified time squence.
#'
#' @description The function computes numerical solutions for a system of ordinary differential equations
#' given an initial state, parameters, and a sequence of time points \eqn{t_1, ..., t_T}.
#'
#' @details This function numerically integrates a system of ODEs using solvers available in the
#' \code{deSolve} package. The ODE system is defined by the user through \code{odefun},
#' which specifies the rate of change for each state variable.
#'
#' The solver computes the values of the state variables at each time step, producing
#' a trajectory of the system's evolution over time.
#'
#' The ODE system must be written in first-order form:
#' \deqn{\frac{dy}{dt} = f(t, y, p)}
#' where:
#' - \eqn{y} is the vector of state variables,
#' - \eqn{t} is time,
#' - \eqn{p} is a set of model parameters,
#' - \eqn{f} is a function returning the derivatives.
#'
#' This function simplifies the process of solving ODEs by managing input formatting
#' and output structure, making it easier to extract and analyze results.
#'
#' For details, see "\href{https://CRAN.R-project.org/package=deSolve}{\code{ode}}".
#'
#' @param odefun function defining the ODE system. It should return a list where the first element is a numeric vector of derivatives.
#' @param times numeric vector of time points where the solution is computed.
#' @param param numeric scalar or vector of parameters to be passed to \code{odefun}.
#' @param init numeric vector or a single-row matrix specifying the initial state values. Each element corresponds to a state variable.
#' @param method character specifying the numerical solver. Default is \code{"lsoda"}.
#' @return A list containing:
#' \itemize{
#'   \item \code{times}: copy of \code{times}
#'   \item \code{<state variable>}: Numeric vectors of computed values for each state variable.
#' }
#'
#' @importFrom deSolve ode
#' @usage dyn_solve(odefun, times, param, init, method)
#' @export
#' @examples
#' library(deSolve)
#' ### Lotka-Volterra equations ###
#' LVmod0D <- function(Time, State, Pars) {
#'   with(as.list(c(State, Pars)),{
#'     IngestC <- rI * P * C
#'     GrowthP <- rG * P * (1 - P/K)
#'     MortC <- rM * C
#'
#'     dP <- GrowthP - IngestC
#'     dC <- IngestC * AE - MortC
#'     return(list(c(dP, dC)))
#'   })
#' }
#'
#' ### Define parameters ###
#' pars <- c(rI = 1.5, rG = 1.5, rM = 2, AE = 1, K = 10)
#'
#' ### Define time sequence ###
#' times <- seq(0, 30, by = 0.01)
#'
#' ### Initial conditions ###
#' state <- c(P = 1, C = 2)
#'
#' ### Solve ODE ###
#' dyn_solve(LVmod0D, times, pars, state)
#'

dyn_solve <- function(odefun, times, param, init, method="lsoda"){
  if (is.null(dim(init))) init <- matrix(init, nrow = 1, dimnames = list(NULL, names(init)))

  d <- ncol(init)
  solv <- apply(init, 1, ode, times = times, func = odefun, parms = param, method=method)[-(1:length(times))]

  dyn_result <- list()
  dyn_result[["times"]] <- times
  for (i in 1:d) {
    dyn_result[[colnames(init)[i]]] <- solv[(length(times)*(i-1)+1):(length(times)*i)]
  }

  return(dyn_result)
}