jacobian {ATNr}

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

O(h^5)

, where

h = eps*bioms

. The choice of eps should ensure that

h^5

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)