`numjac`

Numerical Jacobian matrix of a function, by finite differences.

Category: Other public functions: funfun

Syntax

[dFdy, fac] = numjac(F, t, y, Fty, thresh, fac, vectorized)

Estimates the matrix of partial derivatives of F(t, y) with respect to y.

Description

numjac perturbs one component of y at a time, re-evaluates F and divides the change by the step. This is what stiff ODE solvers need when the user does not supply an analytic Jacobian. The step is relative, so fac should be small, around 1e-7, for an accurate result; a large value gives a coarse slope, not the derivative.

For a linear F the finite difference is exact, up to rounding.

Input arguments

F: Function handle F(t, y) returning a column vector.
t, y: Point where the Jacobian is wanted: a scalar time and a column vector.
Fty: The value F(t, y), already computed, to save one evaluation.
thresh: Vector of lower limits for the size of each component of y when choosing the step.
fac: Vector of relative step sizes; each step is fac(j) times the larger of abs(y(j)) and thresh(j).
vectorized: Nonzero if F can take several columns of y at once.

Output arguments

dFdy: The n-by-n Jacobian, entry (i, j) being the derivative of component i of F with respect to y(j).
fac: Step sizes to reuse on the next call.

Example

Jacobian of a two-variable system

F = @(t, y) [y(1)*y(2); y(1) + y(2)];
y0 = [2; 3];
J = numjac(F, 0, y0, F(0, y0), [1e-6; 1e-6], 1e-7*ones(2, 1), 0);
disp(round(J, 3))

Output

     3     2
     1     1

The exact Jacobian at (2, 3) is [3 2; 1 1].

numjac