ode.R
In targeted: Targeted Inference

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
 collapse = TRUE,
 #dev="png",
 comment = "#>"
)
library(targeted)

## -----------------------------------------------------------------------------
args(targeted::specify_ode)

## -----------------------------------------------------------------------------
dy <- targeted::specify_ode("dy = y - 1;")

## -----------------------------------------------------------------------------
args(targeted::solve_ode)

## -----------------------------------------------------------------------------
t <- seq(0, 10, length.out=1e4)
y <- targeted::solve_ode(dy, t, init=0)
plot(t, y, type='l', lwd=3)

## -----------------------------------------------------------------------------
library(targeted)
ode <- 'dy(0) = p(0)*(y(1)-y(0));
        dy(1) = y(0)*(p(1)-y(2));
        dy(2) = y(0)*y(1)-p(2)*y(2);'
f <- specify_ode(ode)

## -----------------------------------------------------------------------------
tt <- seq(0, 100, length.out=2e4)
y <- solve_ode(f, input=tt, init=c(1, 1, 1), par=c(10, 28, 8/3))
head(y)

## -----------------------------------------------------------------------------
colnames(y) <- c("x","y","z")
scatterplot3d::scatterplot3d(y, cex.symbols=0.1, type='b',
                             color=viridisLite::viridis(nrow(y)))

## -----------------------------------------------------------------------------
n <- 1e4
tt <- seq(0, 10, length.out=n)  # Time
xx <- rep(0, n); xx[(n/3):(2*n/3)] <- 1  # Exogenous input, x(t)
input <- cbind(tt, xx)

## -----------------------------------------------------------------------------
mod <- 'double t = x(0);
        dy = p(0) + p(1)*y + p(2)*x(1)*y + p(3)*x(1)*t;'
dy <- specify_ode(mod)

## -----------------------------------------------------------------------------
yy <- solve_ode(dy, input=input, init=100, c(0, .4, -.5, -5))
plot(tt, yy, type='l', lwd=3, xlab='time', ylab='y')

## -----------------------------------------------------------------------------
sessionInfo()