jacobian function

Estimate the Jacobian matrix of a ODE system

jacobian(bioms, ODE, eps = 1e-06)

Arguments

• bioms: float vector, biomass of species.
• ODE: function that computes the ODEs from one of the model available
• eps: float, scale precision of the numerical approximation.

Returns

A matrix corresponding to the Jacobian of the system estimated at the parameter biomasses

Details

The function provides a numerical estimation of the Jacobian matrix based on the 5 points stencil method. The precision of the method is in

, where

. The choice of eps should ensure that

is always lower to the extinction threshold.

The dimension of the Jacobian matrix are not always matching the number of species in the system. This is because we considered that a perturbation can not correspond to the recolonisation of an extinct species. Therefore, extinct species are removed from the system to calculate the Jacobian matrix.

Examples

library(ATNr)
set.seed(123)
# first run a model to reach equilibrium
masses <- runif(20, 10, 100) #body mass of species
L <- create_Lmatrix(masses, 10, Ropt = 10)
L[L > 0] <- 1
mod <- create_model_Unscaled_nuts(20, 10, 3, masses, L)
mod <- initialise_default_Unscaled_nuts(mod, L)
biomasses <- masses ^ -0.75 * 10 ^ 4 #biomasses of species
biomasses <- append(runif(3, 20, 30), biomasses)
times <- seq(0, 100, 1)
sol <- lsoda_wrapper(times, biomasses, mod)
# get the final biomasses
final.bioms = sol[nrow(sol), -1]
# estimate jacobian
jacobian(final.bioms, mod$ODE)