symengine (version 0.2.6)

dxdt: Solve System of Ordinary Differential Equations

Description

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.

Usage

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)

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.

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`.

Examples

# A differential equation specified with dxdt and ==
x <- Symbol("x")
eq <- dxdt(x) == 1/exp(x)
print(eq)
if (FALSE) {
## Lorenz system
use_vars(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")
}