cvode: cvode

Description Usage Arguments Examples

View source: R/RcppExports.R

Description

CVODE solver to solve stiff ODEs

Usage

1
cvode(time_vector, IC, input_function, reltolerance, abstolerance)

Arguments

time_vector

time vector

IC

Initial Conditions

input_function

External pointer to RHS function

reltolerance

Relative Tolerance (a scalar)

abstolerance

Absolute Tolerance (a vector with length equal to ydot)

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# ODEs described by an R function
ODE_R <- function(t, y){

  # vector containing the right hand side gradients
  ydot = vector(mode = "numeric", length = length(y))

  # R indices start from 1
  ydot[1] = -0.04 * y[1] + 1e04 * y[2] * y[3]
  ydot[3] = 3e07 * 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){

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

                ydot[0] = -0.04 * y[0] + 1e04 * y[1] * y[2];
                ydot[2] = 3e07 * y[1] * y[1];
                ydot[1] = -ydot[0] - ydot[2];

                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)
reltol <- 1e-04
abstol <- c(1e-8,1e-14,1e-6)

## Solving the ODEs using cvode function
df1 <- cvode(time_vec, IC, ODE_R , reltol, abstol)           ## using R
df2 <- cvode(time_vec, IC, ODE_Rcpp , reltol, abstol)        ## using Rcpp

## Check that both solutions are identical
# identical(df1, df2)

sundialr documentation built on July 23, 2018, 1:01 a.m.