numsolve: Numerical solver for Ordinary Differential Equation (ODE) systems.

Description

Calculates numerical solution to the ODE system specified in the ode object and returns the solution path (and optionally its derivatives).

Usage

numsolve(o, time, x0, param, approx_sensitivity = NULL, ...)

Arguments

An object of class ode (created via mak or plk).

time

A numeric vector holding desired time-points at which to evaluate the solution path. By convention of this package, whenever time decreases the parameters and initial conditions are reset, see details.

A non-negative numeric vector or matrix holding the initial conditons for the ODE system.

param

A non-negative numeric vector or matrix holding the parameter values for the ODE system (or a list of these if the ode system has multiple arguments, like ratmak).

approx_sensitivity

Logical or NULL. If NULL (default) the sensitivity equations are not solved and no derivatives are returned. Otherwise if logical, it will solve the sensitivity equations. If TRUE an approximate solution is returned, else a more exact solution (based on the numerical solver) is returned.

...

Additional arguments passed to numeric solver.

Value

If approx_sensitivity is NULL a matrix with the first column holding the time vector and all consecutive columns holding the states. A 0/1 convergence code of the solver is given as an attribute.

If approx_sensitivity is logical a list is returned with

trajectory

A matrix containing the trajectory of the system (as with approx_sensitivity = NULL).

sensitivity_param

A list with arrays of sensitivity solutions of parameters (row = time, colum = state, slice = parameter coordinate).

sensitivity_x0

An array with sensitivity solutions of initial state (row = time, colum = state, slice = initial state coordinate).

For a convergence code of 1, try decreasing tol or increasing step_max in solver object in ode object. Note that the sensitivity equations may have non-finites, even if the convergence code is 0. This is typically due to zero initial states at which the state derivative of the field may be undefined.

Details

As mentioned above, whenever time decreases the system restarts at a new initial condition and with a new parameter vector, called a context/environment. The number of contexts is s, i.e., the number of decreases + 1. To support these s contexts, the new initial conditions and parameter vectors can be supplied through param and x0. For both of these, either a matrix form holding the new initial conditions or parameters in the columns or a concatinated vector will work. In either case param and x0 are coarced to matrices with p (= number of parameters) and d (= dimension of state space) rows respectively (with p and d extracted from o). Hence a check of whether param and x0 have length divisible by p and d respectively is conducted. Both parameter vectors and initial conditions will be recycled if s exceeds the number of these divisors. Therefore, if the same parameter vector and/or initial condition is desired across resets, supply only that vector. A warning will be thrown if p or d is neither a multiple nor a sub-multiple of s.

Examples

Run this code

# NOT RUN {
# Example: Michaelis-Menten system
A <- matrix(
c(1, 1, 0, 0,
  0, 0, 1, 0,
  0, 0, 1, 0), ncol = 4, byrow = TRUE)
B <- matrix(
c(0, 0, 1, 0,
  1, 1, 0, 0,
  1, 0, 0, 1), ncol = 4, byrow = TRUE)
k <- c(1, 2, 0.5)
x0 <- c(E = 1, S = 4, ES = 0, P = 0)
Time <- seq(0, .5, by = .1)

m <- mak(A, B)

# Solution for one context
numsolve(m, Time, x0, k)

# Solution for two contexts (the latter with faster rate)
numsolve(m, c(Time, Time), x0, cbind(k, k * 1.5))

# Solution for two contexts (the latter with different initial condition)
numsolve(m, c(Time, Time), cbind(x0, x0 + 1.5), k)

# As above, but with sensitivity equations are solved (using approximate solution)
numsolve(m, c(Time, Time), cbind(x0, x0 + 1.5), k, TRUE)


# Example: Power law kinetics
A <- matrix(c(1, 0, 1,
              1, 1, 0), byrow = TRUE, nrow = 2)
p <- plk(A)
x0 <- c(10, 4, 1)
theta <- matrix(c(0, -0.25,
                  0.75, 0,
                  0, -0.1), byrow = TRUE, nrow = 3)
numsolve(p, Time, x0, theta)

# }