dxdt {symengine}

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 signature 'DxdtOdeConstructor,ANY'
e1 == e2

ODESystem(
  odesys,
  ...,
  method = "rk5_i",
  atol = 1e-06,
  rtol = 1e-06,
  compile = TRUE
)

## S4 method for signature '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)
## Not run: 
## 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")

## End(Not run)