,
Lawrence Livermore National Security [cph],
Southern Methodist University [cph]| Title: | An Interface to 'SUNDIALS' Ordinary Differential Equation (ODE) Solvers |
|---|---|
| Description: | Provides a way to call the functions in 'SUNDIALS' C ODE solving library (<https://computing.llnl.gov/projects/sundials>). Currently the serial version of ODE solver, 'CVODE', sensitivity calculator 'CVODES' and differential algebraic solver 'IDA' from the 'SUNDIALS' library are implemented. The package requires ODE to be written as an 'R' or 'Rcpp' function and does not require the 'SUNDIALS' library to be installed on the local machine. |
| Authors: | Satyaprakash Nayak [aut, cre, cph]
,
Lawrence Livermore National Security [cph],
Southern Methodist University [cph] |
| Maintainer: | Satyaprakash Nayak <[email protected]> |
| License: | BSD_3_clause + file LICENSE |
| Version: | 0.1.6.2 |
| Built: | 2025-03-15 05:13:43 UTC |
| Source: | https://github.com/sn248/sundialr |
CVODE solver to solve stiff ODEs
cvode( time_vector, IC, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04 )cvode( time_vector, IC, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04 )
time_vector |
time vector |
IC |
Initial Conditions |
input_function |
Right Hand Side function of ODEs |
Parameters |
Parameters input to ODEs |
reltolerance |
Relative Tolerance (a scalar, default value = 1e-04) |
abstolerance |
Absolute Tolerance (a scalar or vector with length equal to ydot (dy/dx), default = 1e-04) |
A data frame. First column is the time-vector, the other columns are values of y in order they are provided.
# Example of solving a set of ODEs with cvode function # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] + p[2]*y[2]*y[3] ydot[2] = p[1]*y[1] - p[2]*y[2]*y[3] - p[3]*y[2]*y[2] ydot[3] = p[3]*y[2]*y[2] # ydot[1] = -0.04 * y[1] + 10000 * y[2] * y[3] # ydot[3] = 30000000 * y[2] * y[2] # ydot[2] = -ydot[1] - ydot[3] ydot } # ODEs can also be described using Rcpp Rcpp::sourceCpp(code = ' #include <Rcpp.h> using namespace Rcpp; // ODE functions defined using Rcpp // [[Rcpp::export]] NumericVector ODE_Rcpp (double t, NumericVector y, NumericVector p){ // Initialize ydot filled with zeros NumericVector ydot(y.length()); ydot[0] = -p[0]*y[0] + p[1]*y[1]*y[2]; ydot[1] = p[0]*y[0] - p[1]*y[1]*y[2] - p[2]*y[1]*y[1]; ydot[2] = p[2]*y[1]*y[1]; return ydot; }') # R code to genrate time vector, IC and solve the equations time_vec <- c(0.0, 0.4, 4.0, 40.0, 4E2, 4E3, 4E4, 4E5, 4E6, 4E7, 4E8, 4E9, 4E10) IC <- c(1,0,0) params <- c(0.04, 10000, 30000000) reltol <- 1e-04 abstol <- c(1e-8,1e-14,1e-6) ## Solving the ODEs using cvode function df1 <- cvode(time_vec, IC, ODE_R , params, reltol, abstol) ## using R df2 <- cvode(time_vec, IC, ODE_Rcpp , params, reltol, abstol) ## using Rcpp ## Check that both solutions are identical # identical(df1, df2)# Example of solving a set of ODEs with cvode function # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] + p[2]*y[2]*y[3] ydot[2] = p[1]*y[1] - p[2]*y[2]*y[3] - p[3]*y[2]*y[2] ydot[3] = p[3]*y[2]*y[2] # ydot[1] = -0.04 * y[1] + 10000 * y[2] * y[3] # ydot[3] = 30000000 * y[2] * y[2] # ydot[2] = -ydot[1] - ydot[3] ydot } # ODEs can also be described using Rcpp Rcpp::sourceCpp(code = ' #include <Rcpp.h> using namespace Rcpp; // ODE functions defined using Rcpp // [[Rcpp::export]] NumericVector ODE_Rcpp (double t, NumericVector y, NumericVector p){ // Initialize ydot filled with zeros NumericVector ydot(y.length()); ydot[0] = -p[0]*y[0] + p[1]*y[1]*y[2]; ydot[1] = p[0]*y[0] - p[1]*y[1]*y[2] - p[2]*y[1]*y[1]; ydot[2] = p[2]*y[1]*y[1]; return ydot; }') # R code to genrate time vector, IC and solve the equations time_vec <- c(0.0, 0.4, 4.0, 40.0, 4E2, 4E3, 4E4, 4E5, 4E6, 4E7, 4E8, 4E9, 4E10) IC <- c(1,0,0) params <- c(0.04, 10000, 30000000) reltol <- 1e-04 abstol <- c(1e-8,1e-14,1e-6) ## Solving the ODEs using cvode function df1 <- cvode(time_vec, IC, ODE_R , params, reltol, abstol) ## using R df2 <- cvode(time_vec, IC, ODE_Rcpp , params, reltol, abstol) ## using Rcpp ## Check that both solutions are identical # identical(df1, df2)
CVODES solver to solve ODEs and calculate sensitivities
cvodes( time_vector, IC, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04, SensType = "STG", ErrCon = "F" )cvodes( time_vector, IC, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04, SensType = "STG", ErrCon = "F" )
time_vector |
time vector |
IC |
Initial Conditions |
input_function |
Right Hand Side function of ODEs |
Parameters |
Parameters input to ODEs |
reltolerance |
Relative Tolerance (a scalar, default value = 1e-04) |
abstolerance |
Absolute Tolerance (a scalar or vector with length equal to ydot, default = 1e-04) |
SensType |
Sensitivity Type - allowed values are "STG" (for Staggered, default) or "SIM" (for Simultaneous) |
ErrCon |
Error Control - allowed values are TRUE or FALSE (default) |
A data frame. First column is the time-vector, the next y * p columns are sensitivities of y1 w.r.t all parameters, then y2 w.r.t all parameters etc. y is the state vector, p is the parameter vector
# Example of solving a set sensitivity equations for ODEs with cvodes function # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] + p[2]*y[2]*y[3] ydot[2] = p[1]*y[1] - p[2]*y[2]*y[3] - p[3]*y[2]*y[2] ydot[3] = p[3]*y[2]*y[2] # ydot[1] = -0.04 * y[1] + 10000 * y[2] * y[3] # ydot[3] = 30000000 * y[2] * y[2] # ydot[2] = -ydot[1] - ydot[3] ydot } # ODEs can also be described using Rcpp Rcpp::sourceCpp(code = ' #include <Rcpp.h> using namespace Rcpp; // ODE functions defined using Rcpp // [[Rcpp::export]] NumericVector ODE_Rcpp (double t, NumericVector y, NumericVector p){ // Initialize ydot filled with zeros NumericVector ydot(y.length()); ydot[0] = -p[0]*y[0] + p[1]*y[1]*y[2]; ydot[1] = p[0]*y[0] - p[1]*y[1]*y[2] - p[2]*y[1]*y[1]; ydot[2] = p[2]*y[1]*y[1]; return ydot; }') # R code to genrate time vector, IC and solve the equations time_vec <- c(0.0, 0.4, 4.0, 40.0, 4E2, 4E3, 4E4, 4E5, 4E6, 4E7, 4E8, 4E9, 4E10) IC <- c(1,0,0) params <- c(0.04, 10000, 30000000) reltol <- 1e-04 abstol <- c(1e-8,1e-14,1e-6) ## Solving the ODEs and Sensitivities using cvodes function df1 <- cvodes(time_vec, IC, ODE_R , params, reltol, abstol,"STG",FALSE) ## using R df2 <- cvodes(time_vec, IC, ODE_Rcpp , params, reltol, abstol,"STG",FALSE) ## using Rcpp ## Check that both solutions are identical # identical(df1, df2)# Example of solving a set sensitivity equations for ODEs with cvodes function # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] + p[2]*y[2]*y[3] ydot[2] = p[1]*y[1] - p[2]*y[2]*y[3] - p[3]*y[2]*y[2] ydot[3] = p[3]*y[2]*y[2] # ydot[1] = -0.04 * y[1] + 10000 * y[2] * y[3] # ydot[3] = 30000000 * y[2] * y[2] # ydot[2] = -ydot[1] - ydot[3] ydot } # ODEs can also be described using Rcpp Rcpp::sourceCpp(code = ' #include <Rcpp.h> using namespace Rcpp; // ODE functions defined using Rcpp // [[Rcpp::export]] NumericVector ODE_Rcpp (double t, NumericVector y, NumericVector p){ // Initialize ydot filled with zeros NumericVector ydot(y.length()); ydot[0] = -p[0]*y[0] + p[1]*y[1]*y[2]; ydot[1] = p[0]*y[0] - p[1]*y[1]*y[2] - p[2]*y[1]*y[1]; ydot[2] = p[2]*y[1]*y[1]; return ydot; }') # R code to genrate time vector, IC and solve the equations time_vec <- c(0.0, 0.4, 4.0, 40.0, 4E2, 4E3, 4E4, 4E5, 4E6, 4E7, 4E8, 4E9, 4E10) IC <- c(1,0,0) params <- c(0.04, 10000, 30000000) reltol <- 1e-04 abstol <- c(1e-8,1e-14,1e-6) ## Solving the ODEs and Sensitivities using cvodes function df1 <- cvodes(time_vec, IC, ODE_R , params, reltol, abstol,"STG",FALSE) ## using R df2 <- cvodes(time_vec, IC, ODE_Rcpp , params, reltol, abstol,"STG",FALSE) ## using Rcpp ## Check that both solutions are identical # identical(df1, df2)
CVSOLVE solver to solve stiff ODEs with discontinuties
cvsolve( time_vector, IC, input_function, Parameters, Events = NULL, reltolerance = 1e-04, abstolerance = 1e-04 )cvsolve( time_vector, IC, input_function, Parameters, Events = NULL, reltolerance = 1e-04, abstolerance = 1e-04 )
time_vector |
time vector |
IC |
Initial Conditions |
input_function |
Right Hand Side function of ODEs |
Parameters |
Parameters input to ODEs |
Events |
Discontinuities in the solution (a DataFrame, default value is NULL) |
reltolerance |
Relative Tolerance (a scalar, default value = 1e-04) |
abstolerance |
Absolute Tolerance (a scalar or vector with length equal to ydot, default = 1e-04) |
A data frame. First column is the time-vector, the other columns are values of y in order they are provided.
# Example of solving a set of ODEs with multiple discontinuities using cvsolve # A simple One dimensional equation, y = -0.1 * y # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] ydot } # R code to generate time vector, IC and solve the equations TSAMP <- seq(from = 0, to = 100, by = 0.1) # sampling time points IC <- c(1) params <- c(0.1) # A dataset describing the dosing at times at which additions to y[1] are to be done # Names of the columns don't matter, but they MUST be in the order of state index, # times and Values at discontinuity. TDOSE <- data.frame(ID = 1, TIMES = c(0, 10, 20, 30, 40, 50), VAL = 100) df1 <- cvsolve(TSAMP, c(1), ODE_R, params) # solving without any discontinuity df2 <- cvsolve(TSAMP, c(1), ODE_R, params, TDOSE) # solving with discontinuity# Example of solving a set of ODEs with multiple discontinuities using cvsolve # A simple One dimensional equation, y = -0.1 * y # ODEs described by an R function ODE_R <- function(t, y, p){ # vector containing the right hand side gradients ydot = vector(mode = "numeric", length = length(y)) # R indices start from 1 ydot[1] = -p[1]*y[1] ydot } # R code to generate time vector, IC and solve the equations TSAMP <- seq(from = 0, to = 100, by = 0.1) # sampling time points IC <- c(1) params <- c(0.1) # A dataset describing the dosing at times at which additions to y[1] are to be done # Names of the columns don't matter, but they MUST be in the order of state index, # times and Values at discontinuity. TDOSE <- data.frame(ID = 1, TIMES = c(0, 10, 20, 30, 40, 50), VAL = 100) df1 <- cvsolve(TSAMP, c(1), ODE_R, params) # solving without any discontinuity df2 <- cvsolve(TSAMP, c(1), ODE_R, params, TDOSE) # solving with discontinuity
IDA solver to solve stiff DAEs
ida( time_vector, IC, IRes, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04 )ida( time_vector, IC, IRes, input_function, Parameters, reltolerance = 1e-04, abstolerance = 1e-04 )
time_vector |
time vector |
IC |
Initial Value of y |
IRes |
Inital Value of ydot |
input_function |
Right Hand Side function of DAEs |
Parameters |
Parameters input to ODEs |
reltolerance |
Relative Tolerance (a scalar, default value = 1e-04) |
abstolerance |
Absolute Tolerance (a scalar or vector with length equal to ydot, default = 1e-04) |
A data frame. First column is the time-vector, the other columns are values of y in order they are provided.