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 nonzero elements are restricted to a number of bands above and below the diagonal.
1 2  jacobian.band(y, func, bandup = 1, banddown = 1,
dy = NULL, time = 0, parms = NULL, pert = 1e8, ...)

y 
(state) variables, a vector; if 
func 
function that calculates one function value for each element
of 
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 
time 
time, passed to function 
parms 
parameter values, passed to function 
pert 
numerical perturbation factor; increase depending on precision of model solution. 
... 
other arguments passed to function 
The function func
that estimates the rate of change of the state
variables has to be consistent with functions called from Rpackage
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 xvalues.
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 
[,4]  0  0  9  10 
the banded jacobian will be:
[,1],  [,2],  [,3],  [,4]  
[,1]  0  2  5  8 
[,2]  1  4  7  10 
[,3]  3  6  9  0 
Karline Soetaert <karline.soetaert@nioz.nl>
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
1 2 3 4 5 6 7 8 9 10 11  ## =======================================================================
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)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.