jacobian {ATNr} | R Documentation |
Estimate the Jacobian matrix of a ODE system
Description
Estimate the Jacobian matrix of a ODE system
Usage
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. |
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.
Value
A matrix corresponding to the Jacobian of the system estimated at the parameter biomasses
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)