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

1
2
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.

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

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

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
[,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

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

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

 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? or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.