This is a wrapper of the odeintr R package using
symengine objects to specify the ODE system and C code
generation functionality from symengine to generate the
C++ source. The dxdt function and defined == S4 method
allow one to intuitively specify the ODE system with symengine
objects. The ODESystem will generate C++ source
and compile on the fly with Rcpp. Then predict can be used to get
results.
dxdt(x)
# S4 method for DxdtOdeConstructor,ANY
==(e1, e2)
ODESystem(
odesys,
...,
method = "rk5_i",
atol = 1e-06,
rtol = 1e-06,
compile = TRUE
)
# S4 method for ODESystem
predict(object, init, duration, step_size = 1, start = 0)Arguments
- x
A SymEngine Basic object of type Symbol or a R object that will be converted to
Symbol(x).- e1
A DxdtOdeConstructor S4 object which can be returned by `dxdt`.
- e2
A Basic object or an R object that will be converted to `S(e2)`.
- odesys, ...
DxdtOde S4 objects that can be returned with `dxdt(x) == rhs`. Or `odesys` can be a list of DxdtOde S4 objects when there is no dot arguments.
- method, atol, rtol
Passed to `odeintr::compile_sys`.
- compile
Logical, whether to compile the C++ source. Useful if you only want to obtain the code.
- object
A ODESystem S4 object.
- init
A numeric vector specifying the initial conditions. It can be named with the variable names or it can be unnamed but in the same of order of equations.
- duration, step_size, start
Passed to the function generated by `odeintr::compile_sys`.
Value
dxdt returns a DxdtOdeConstructor S4 object.
S4 method of `==` for "DxdtOdeConstructor" returns a DxdtOde S4 object.
`ODESystem` returns a ODESystem S4 object.
`predict` returns a dataframe.
Examples
# A differential equation specified with dxdt and ==
x <- Symbol("x")
eq <- dxdt(x) == 1/exp(x)
print(eq)
#> Ordinary differential equation:
#> d(x)/dt == 1.0*exp(-x)
# \donttest{
## Lorenz system
use_vars(x, y, z)
#> Initializing ‘x’, ‘y’, ‘z’
sigma <- 10
rho <- 28
beta <- 8/3
lorenz_sys <- ODESystem(
dxdt(x) == sigma * (y - x),
dxdt(y) == (rho - z) * x - y,
dxdt(z) == - beta * z + x * y
)
res <- predict(
lorenz_sys, init = c(x = 1, y = 1, z = 1), duration = 100, step_size = 0.001
)
plot(res[, c(2, 4)], type = 'l', col = "steelblue", main = "Lorenz Attractor")
# }