ida

IDA solver to solve stiff DAEs

Usage

ida(
  time_vector,
  IC,
  IRes,
  input_function,
  Parameters,
  reltolerance = 1e-04,
  abstolerance = 1e-04,
  jacobian = NULL
)

Arguments

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)

jacobian

(Optional) Jacobian with signature function(t, y, ydot, cj, p) returning an n-by-n matrix of dF/dy + cj*dF/dydot. Default NULL.

Value

A Matrix. First column is the time-vector, the other columns are values of y in order they are provided.

Examples

# Example of solving a set of Differential Algebraic Equations (DAEs)
# with IDA function
# DAEs (residuals) described by an R function
DAE_R <- function(t, y, ydot, p){

  # vector containing the residuals
  res = vector(mode = "numeric", length = length(y))

  # R indices start from 1
  res[1] <- -0.04 * y[1] + 10000 * y[2] * y[3] - ydot[1]
  res[2] <- -res[1] - 30000000 * y[2] * y[2] - ydot[2]
  res[3] <- y[1] + y[2] + y[3] - 1.0

  res
}

# DAEs can also be described using Rcpp
Rcpp::sourceCpp(code = '

                #include <Rcpp.h>
                using namespace Rcpp;

                // ODE functions defined using Rcpp
                // [[Rcpp::export]]
                NumericVector DAE_Rcpp (double t, NumericVector y, NumericVector ydot, NumericVector p){

                // Initialize ydot filled with zeros
                NumericVector res(y.length());

                res[0] = -0.04 * y[0] + 10000 * y[1] * y[2];
                res[1] = -res[0] - 30000000 * y[1] * y[1] - ydot[1];
                res[0] = res[0] - ydot[0];
                res[2] = y[0] + y[1] + y[2] - 1.0;

                return res;

                }')

# 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)
IRes <- c(-0.4, 0.4, 0)
params <- c(0.04, 10000, 30000000)
reltol <- 1e-04
abstol <- c(1e-8,1e-14,1e-6)

## Solving the DAEs using the ida function
df1 <- ida(time_vec, IC, IRes, DAE_R , params, reltol, abstol)           ## using R
df2 <- ida(time_vec, IC, IRes, DAE_Rcpp , params, reltol, abstol)        ## using Rcpp

## Solving with a manual Jacobian
## J[i,j] = dF_i/dy_j + cj * dF_i/dydot_j
##
## F1 = -0.04*y1 + 1e4*y2*y3 - ydot1
## F2 = 0.04*y1 - 1e4*y2*y3 - 3e7*y2^2 - ydot2
## F3 = y1 + y2 + y3 - 1  (algebraic constraint)
DAE_jac <- function(t, y, ydot, p) {
  res <- numeric(length(y))
  f1     <- -0.04 * y[1] + 10000 * y[2] * y[3]
  res[1] <- f1 - ydot[1]
  res[2] <- -f1 - 30000000 * y[2] * y[2] - ydot[2]
  res[3] <- y[1] + y[2] + y[3] - 1.0
  res
}
JAC_IDA <- function(t, y, ydot, cj, p) {
  matrix(c(
    -0.04 - cj,   0.04,    1,
    10000*y[3],  -10000*y[3] - 60000000*y[2] - cj,   1,
    10000*y[2],  -10000*y[2],                         1
  ), nrow = 3, ncol = 3)
}
df3 <- ida(time_vec, IC, IRes, DAE_jac, params, reltol, abstol, jacobian = JAC_IDA)