# numsolve: Numerical solver for Ordinary Differential Equation (ODE)... In episode: Estimation with Penalisation in Systems of Ordinary Differential Equations

## Description

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

## Usage

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

## Arguments

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

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

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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37``` ```# 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) ```