Ordinary differential equation solver using lsoda (C++ code)

Usage

ode_cpp(y, times, func, rtol = 1e-06, atol = 1e-06)

Arguments

y: vector of initial state values
times: vector of times – including the start time
func: R function with signature function(t,y) that returns a list: the first list element is a vector for dy/dt; the second list element, if it exists, is a vector of result calculations to be retained.
rtol: double for the relative tolerance
atol: double for the absolute tolerance

Value

a matrix for times in the first column and the state andd results values in the other columns.

Examples

  times = c(0,0.4*10^(0:10))
 y = c(1,0,0)
 func = function(t,y) {
     ydot = rep(0,3)
     ydot[1] = 1.0E4 * y[2] * y[3] - .04E0 * y[1]
     ydot[3] = 3.0E7 * y[2] * y[2]
     ydot[2] = -1.0 * (ydot[1] + ydot[3])
     list(ydot, sum(y))
 }
 lsoda::ode_cpp(y,times,func, rtol=1e-8, atol=1e-8)
#>        time           y1           y2         y3 res1
#>  [1,] 0e+00 1.000000e+00 0.000000e+00 0.00000000    1
#>  [2,] 4e-01 9.851721e-01 3.386395e-05 0.01479403    1
#>  [3,] 4e+00 9.055187e-01 2.240475e-05 0.09445893    1
#>  [4,] 4e+01 7.158270e-01 9.185532e-06 0.28416381    1
#>  [5,] 4e+02 4.505186e-01 3.222901e-06 0.54947815    1
#>  [6,] 4e+03 1.832023e-01 8.942356e-07 0.81679679    1
#>  [7,] 4e+04 3.898338e-02 1.621768e-07 0.96101646    1
#>  [8,] 4e+05 4.938305e-03 1.985006e-08 0.99506168    1
#>  [9,] 4e+06 5.168319e-04 2.068387e-09 0.99948317    1
#> [10,] 4e+07 5.204267e-05 2.081814e-10 0.99994796    1
#> [11,] 4e+08 5.214038e-06 2.085626e-11 0.99999479    1
#> [12,] 4e+09 5.280602e-07 2.112242e-12 0.99999947    1