rootSolve (version 1.8.2.1)

jacobian.band: Banded jacobian matrix for a system of ODEs (ordinary differential equations)

Description

Given a vector of (state) variables y, and a function that estimates a function value for each (state) variable (e.g. the rate of change), estimates the Jacobian matrix (\(d(f(y))/d(y)\)).

Assumes a banded structure of the Jacobian matrix, i.e. where the non-zero elements are restricted to a number of bands above and below the diagonal.

Usage

jacobian.band(y, func, bandup = 1, banddown = 1, 
              dy = NULL, time = 0, parms = NULL, pert = 1e-8, ...)

Arguments

y

(state) variables, a vector; if y has a name attribute, the names will be used to label the jacobian matrix columns.

func

function that calculates one function value for each element of y; if an ODE system, func calculates the rate of change (see details).

bandup

number of nonzero bands above the diagonal of the Jacobian matrix.

banddown

number of nonzero bands below the diagonal of the Jacobian matrix.

dy

reference function value; if not specified, it will be estimated by calling func.

time

time, passed to function func.

parms

parameter values, passed to function func.

pert

numerical perturbation factor; increase depending on precision of model solution.

...

other arguments passed to function func.

Value

Jacobian matrix, in banded format, i.e. only the nonzero bands near the diagonal form the rows of the Jacobian.

this matrix has bandup+banddown+1 rows, while the number of columns equal the length of y.

Thus, if the full Jacobian is given by:

	[,1],	[,2],	[,3],	[,4]
[,1]	1	2	0	0
[,2]	3	4	5	0
[,3]	0	6	7	8

the banded jacobian will be:

	[,1],	[,2],	[,3],	[,4]
[,1]	0	2	5	8
[,2]	1	4	7	10

Details

The function func that estimates the rate of change of the state variables has to be consistent with functions called from R-package deSolve, which contains integration routines.

This function call is as: function(time,y,parms,...) where

• y : (state) variable values at which the Jacobian is estimated.

• parms: parameter vector - need not be used.

• time: time at which the Jacobian is estimated - in general, time will not be used.

• ...: (optional) any other arguments

The Jacobian is estimated numerically, by perturbing the x-values.

See Also

jacobian.full, estimates the Jacobian matrix assuming a full matrix.

hessian, estimates the Hessian matrix.

gradient, for a full (not necessarily square) gradient matrix and where the function call is simpler.

uniroot.all, to solve for all roots of one (nonlinear) equation

multiroot, to solve n roots of n (nonlinear) equations

Examples

# NOT RUN {
## =======================================================================

mod <- function (t = 0, y, parms = NULL,...) {
  dy1 <-  y[1] + 2*y[2]
  dy2 <-3*y[1] + 4*y[2] + 5*y[3]
  dy3 <-         6*y[2] + 7*y[3] + 8*y[4]
  dy4 <-                  9*y[3] +10*y[4]
  return(as.list(c(dy1, dy2, dy3, dy4)))
}

jacobian.band(y = c(1, 2, 3, 4), func = mod)
# }