In this vignette we explain the numerical scheme used in mizer. We will not go into the details of the model itself, which is described in the model description vignette. We will focus on how the model is discretised and how the resulting difference equations are solved.
The weight grid
The model dynamics are described by partial differential equations (PDEs) for the species number densities \(N_i(w)\) and the resource number density \(N_R(w)\). These depend on the individual weight \(w\). To solve these equations numerically, we discretise the weight axis.
We choose a grid of weights \(w_1, w_2, \ldots, w_K\). The grid is logarithmically spaced, meaning that the ratio of consecutive weights is constant: \[ w_{j+1} / w_j = 10^{\Delta x} = \text{const}. \] The logarithmic spacing is chosen because the weights of fish span many orders of magnitude, from milligrams to megagrams.
This grid defines a set of bins. The \(j\)-th bin is the interval \([w_j, w_{j+1})\). The width of the \(j\)-th bin is \[ \Delta w_j = w_{j+1} - w_j = w_j (10^{\Delta x} - 1). \] Note that \(\Delta w_j\) increases with \(w_j\), which is appropriate for a logarithmic grid.
The number density \(N_i(w)\) is approximated by a constant value \(N_{i,j}\) within each bin. Thus \(N_{i,j}\) represents the density at the start of the bin \(w_j\). The number of individuals in the \(j\)-th bin is \(N_{i,j} \Delta w_j\).
The Transport Equation
The time evolution of the number density \(N_i(w)\) is described by the McKendrick-von Foerster equation with an added diffusion term: \[ \frac{\partial N_i}{\partial t} + \frac{\partial}{\partial w} \left( g_i N_i - \frac{1}{2}\frac{\partial(d_i N_i)}{\partial w} \right) = -\mu_i N_i \] where \(g_i(w)\) is the somatic growth rate, \(d_i(w)\) is the diffusion coefficient and \(\mu_i(w)\) is the total mortality rate.
We discretise this equation using a finite volume scheme.
Discretisation of the Fluxes
The term inside the derivative with respect to \(w\) is the flux \(J_i(w)\): \[ J_i(w) = g_i(w) N_i(w) - \frac{1}{2}\frac{\partial(d_i(w) N_i(w))}{\partial w} \] We consider the \(j\)-th size bin \([w_j, w_{j+1}]\). Integrating the conservation equation over this bin gives: \[ \Delta w_j \frac{\partial N_{i,j}}{\partial t} + J_i(w_{j+1}) - J_i(w_j) = -\int_{w_j}^{w_{j+1}} \mu_i(w) N_i(w) dw \] Approximating the integral and dividing by \(\Delta w_j\): \[ \frac{N_{i,j}^{t+1} - N_{i,j}^t}{\Delta t} + \frac{J_{i,j+1} - J_{i,j}}{\Delta w_j} = -\mu_i(w_j) N_{i,j}^{t+1} \] where \(J_{i,j}\) represents the flux at the boundary \(w_j\).
Advective Flux
For the advective part of the flux \(g_i N_i\), we use an upwind approximation. Since fish grow from smaller to larger sizes (\(g_i > 0\)), the flux at boundary \(w_j\) is determined by the density in the bin below (\(j-1\)): \[ J_{i,j}^{adv} = g_i(w_{j-1}) N_{i,j-1} \] (Note: for \(j=j_{min}\), this boundary flux is the recruitment \(R_{dd,i}\)).
Diffusive Flux
For the diffusive part of the flux \(-\frac{1}{2}\frac{\partial(d_i N_i)}{\partial w}\), we use a central difference approximation at the boundary \(w_j\). We approximate the gradient using the densities in the adjacent bins \(j-1\) and \(j\): \[ J_{i,j}^{diff} \approx -\frac{1}{2} \frac{(d_i N_i)_j - (d_i N_i)_{j-1}}{w_j - w_{j-1}} = -\frac{1}{2} \frac{d_i(w_j) N_{i,j} - d_i(w_{j-1}) N_{i,j-1}}{\Delta w_{j-1}} \] Note the use of \(N_{i,j}\) (density in bin \(j\)) and \(N_{i,j-1}\) (density in bin \(j-1\)) to estimate the value at the interface \(w_j\).
The total flux at boundary \(w_j\) is: \[ J_{i,j} = g_i(w_{j-1}) N_{i,j-1} - \frac{1}{2} \frac{d_i(w_j) N_{i,j} - d_i(w_{j-1}) N_{i,j-1}}{\Delta w_{j-1}} \] And similarly at boundary \(w_{j+1}\): \[ J_{i,j+1} = g_i(w_j) N_{i,j} - \frac{1}{2} \frac{d_i(w_{j+1}) N_{i,j+1} - d_i(w_j) N_{i,j}}{\Delta w_j} \]
Limitation on Time Step
With the diffusion term, an explicit time discretisation would require a very small time step for stability (\(\Delta t \sim \Delta w^2\)). Therefore, we use a semi-implicit scheme where the densities \(N_i\) are evaluated at time \(t+1\), but the rates (\(g_i, \mu_i, d_i\)) are evaluated at time \(t\).
Discretised Equation
Substituting the fluxes into the conservation equation: \[ \frac{N_{i,j}^{t+1} - N_{i,j}^t}{\Delta t} + \frac{1}{\Delta w_j} \left( J_{i,j+1}^{t+1} - J_{i,j}^{t+1} \right) = -\mu_i(w_j) N_{i,j}^{t+1} \] This leads to a linear system of the form: \[ A_{i,j} N_{i,j-1}^{t+1} + B_{i,j} N_{i,j}^{t+1} + C_{i,j} N_{i,j+1}^{t+1} = S_{i,j} \] where \(S_{i,j} = N_{i,j}^t\). This is a tridiagonal system for each species \(i\), which can be solved efficiently (e.g., using the Thomas algorithm).
The coefficients are: \[ \begin{aligned} A_{i,j} &= -\frac{\Delta t}{\Delta w_j} \left( g_i(w_{j-1}) + \frac{1}{2} \frac{d_i(w_{j-1})}{\Delta w_{j-1}} \right) \\ C_{i,j} &= -\frac{\Delta t}{\Delta w_j} \left( \frac{1}{2} \frac{d_i(w_{j+1})}{\Delta w_j} \right) \\ B_{i,j} &= 1 + \Delta t \mu_i(w_j) + \frac{\Delta t}{\Delta w_j} \left( g_i(w_j) + \frac{1}{2} \frac{d_i(w_j)}{\Delta w_j} + \frac{1}{2} \frac{d_i(w_j)}{\Delta w_{j-1}} \right) \end{aligned} \]
Boundary Conditions
At the smallest size (\(j=j_{min}\)): The flux entering the grid is determined by recruitment. \[ J_{i, j_{min}} = R_{dd, i} \] (We assume diffusive flux at the lower boundary is negligible or incorporated into \(R_{dd}\)). The equation for the first bin becomes: \[ \frac{N_{i,j_{min}}^{t+1} - N_{i,j_{min}}^t}{\Delta t} + \frac{J_{i, j_{min}+1}^{t+1} - R_{dd,i}}{\Delta w_{j_{min}}} = -\mu_i(w_{j_{min}}) N_{i,j_{min}}^{t+1} \] This involves \(N_{i, j_{min}}^{t+1}\) and \(N_{i, j_{min}+1}^{t+1}\). Comparing this to the general discretised equation translates to modifying the first row (\(j=j_{min}\)) of our tri-diagonal matrices:
- The \(j_{min}-1\) term does not exist, so \(A_{i,j_{min}} = 0\).
- The upward diffusion term from below the boundary is omitted, so \(B_{i, j_{min}}\) does not have the \(\frac{1}{2} \frac{d_i(w_{j_{min}})}{\Delta w_{j_{min}-1}}\) component: \[ B_{i,j_{min}} = 1 + \Delta t \mu_i(w_{j_{min}}) + \frac{\Delta t}{\Delta w_{j_{min}}} \left( g_i(w_{j_{min}}) + \frac{1}{2} \frac{d_i(w_{j_{min}})}{\Delta w_{j_{min}}} \right) \]
- The coefficient \(C_{i, j_{min}}\) remains unchanged from the general formula.
- The recruitment flux enters as a source term, so it is added to the right-hand side \(S_{i, j_{min}}\): \[ S_{i,j_{min}} = N_{i,j_{min}}^t + \frac{\Delta t}{\Delta w_{j_{min}}} R_{dd, i} \]
(Additionally, for any size classes below the recruitment size \(j < j_{min}\), we set all coefficients in the matrices \(A\), \(B\), \(C\) and vector \(S\) to \(0\) to avoid any dynamics in that range).
At the largest size (\(j=j_{max}\)): We typically assume that densities drop to zero beyond the maximum size, \(N_{i, j_{max}+1} = 0\). The flux leaving the grid is: \[ J_{i, j_{max}+1} = g_i(w_{j_{max}}) N_{i, j_{max}} - \frac{1}{2} \frac{0 - d_i(w_{j_{max}}) N_{i, j_{max}}}{\Delta w_{j_{max}}} \] This means the term \(C_{i, j_{max}}\) multiplying \(N_{i, j_{max}+1}^{t+1}\) is not needed, so \(C_{i, j_{max}} = 0\). The coefficients \(A_{i, j_{max}}\) and \(B_{i, j_{max}}\) use the standard formulas.
Numerical Diffusion
The upwind scheme used for the advective term introduces numerical diffusion. This is a well-known property of first-order upwind schemes. We can estimate the magnitude of this diffusion by expanding the discretised term using a Taylor series.
The discretised equation for the transport (advection only, with constant rates for simplicity) is: \[ \frac{N_j^{t+1} - N_j^t}{\Delta t} + g \frac{N_j^{t+1} - N_{j-1}^{t+1}}{\Delta w} = 0 \] Expanding \(N(w, t)\) around \((w_j, t+\Delta t)\) leads to the following leading order error terms: \[ \frac{\partial N}{\partial t} + g \frac{\partial N}{\partial w} = \frac{g \Delta w}{2} \left( 1 + \frac{g \Delta t}{\Delta w} \right) \frac{\partial^2 N}{\partial w^2} \] The coefficient of the second derivative represents the numerical diffusivity: \[ D_{num} = \frac{g \Delta w}{2} (1 + C) \] where \(C = \frac{g \Delta t}{\Delta w}\) is the Courant-Friedrichs-Lewy (CFL) number. Comparing this to the Mizer diffusion equation form (where the diffusion term is \(\frac{\partial}{\partial w} ( \frac{1}{2} \frac{\partial (D N)}{\partial w} )\)), the effective diffusion parameter is: \[ d_{num}(w) \approx g(w) \Delta w (1 + C(w)) \] Since \(\Delta w \approx w \ln(\beta)\), this is: \[ d_{num}(w) \approx g(w) w \ln(\beta) \left( 1 + \frac{g(w) \Delta t}{w \ln(\beta)} \right) = g(w) w \ln(\beta) + g(w)^2 \Delta t \] This means the numerical scheme behaves as if there is a diffusion \(d_{num}\). This numerical diffusion has two components: one from spatial discretisation (scaling with \(\Delta w\)) and one from time stepping (scaling with \(\Delta t\)).
Steady-State Solution
When solving the steady-state ODE instead of the time-dependent PDE, we are looking for a state where the population densities do not change over time, meaning \(N_{i,j}^{t+1} = N_{i,j}^t = N_{i,j}^*\).
Substituting this into our discretised linear system: \[ A_{i,j} N_{i,j-1}^* + B_{i,j} N_{i,j}^* + C_{i,j} N_{i,j+1}^* = S_{i,j} \] Recall that for \(j > j_{min}\), \(S_{i,j} = N_{i,j}^t\). The equation simplifies to: \[ A_{i,j} N_{i,j-1}^* + (B_{i,j} - 1) N_{i,j}^* + C_{i,j} N_{i,j+1}^* = 0 \]
To find the steady-state population densities \(N^*\), we formulate a new time-independent tridiagonal system: \[ \tilde{A}_{i,j} N_{i,j-1}^* + \tilde{B}_{i,j} N_{i,j}^* + \tilde{C}_{i,j} N_{i,j+1}^* = \tilde{S}_{i,j} \] To eliminate the explicit dependence on the time step \(\Delta t\), we can divide the equation by \(\Delta t\). The modified coefficients \(\tilde{A}, \tilde{B}, \tilde{C}\) defining the new tri-diagonal system are: \[ \begin{aligned} \tilde{A}_{i,j} &= \frac{A_{i,j}}{\Delta t} = -\frac{1}{\Delta w_j} \left( g_i(w_{j-1}) + \frac{1}{2} \frac{d_i(w_{j-1})}{\Delta w_{j-1}} \right) \\ \tilde{C}_{i,j} &= \frac{C_{i,j}}{\Delta t} = -\frac{1}{\Delta w_j} \left( \frac{1}{2} \frac{d_i(w_{j+1})}{\Delta w_j} \right) \\ \tilde{B}_{i,j} &= \frac{B_{i,j} - 1}{\Delta t} = \mu_i(w_j) + \frac{1}{\Delta w_j} \left( g_i(w_j) + \frac{1}{2} \frac{d_i(w_j)}{\Delta w_j} + \frac{1}{2} \frac{d_i(w_j)}{\Delta w_{j-1}} \right) \end{aligned} \] Notice that \(\tilde{A}_{i,j}\) and \(\tilde{C}_{i,j}\) are exactly the expressions for \(A_{i,j}\) and \(C_{i,j}\) evaluated at \(\Delta t = 1\). Similarly, \(\tilde{B}_{i,j}\) is exactly the expression for \(B_{i,j} - 1\) evaluated at \(\Delta t = 1\).
Boundary conditions for the steady state:
For the smallest size (\(j=j_{min}\)), the original equation had a source term due to recruitment: \[ S_{i, j_{min}} = N_{i,j_{min}}^t + \frac{\Delta t}{\Delta w_{j_{min}}} R_{dd, i} \] Following the same logic of setting \(N^{t+1} = N^t = N^*\) and dividing by \(\Delta t\), the right-hand side vector \(\tilde{S}_{i,j}\) for the steady-state system becomes purely the recruitment flux term. If we again observe the original term \(\frac{\Delta t}{\Delta w_{j_{min}}} R_{dd, i}\) when evaluated at \(\Delta t = 1\), we get our new source vector: \[ \tilde{S}_{i, j_{min}} = \frac{R_{dd, i}}{\Delta w_{j_{min}}} \] For all other \(j > j_{min}\), \(\tilde{S}_{i,j} = 0\).
The boundary condition modifications at the edges of the grid remain the same conceptually: \(\tilde{A}_{i,j_{min}} = 0\), the upward diffusion term is omitted from \(\tilde{B}_{i, j_{min}}\), and \(\tilde{C}_{i, j_{max}} = 0\). For any \(j < j_{min}\), all matrix entries remain zero.
In code, this means that the steady-state coefficients for the matrix multiplication (\(\tilde{A}, \tilde{B}, \tilde{C}\)) and the constant vector (\(\tilde{S}\)) can be calculated by calling the standard coefficient function but simply setting \(\Delta t = 1\), and dropping the \(+1\) and \(+N_{i,j}^t\) from the resulting \(B\) and \(S\) variables respectively.
With these modified matrices, the steady-state densities can be calculated directly by solving the linear system avoiding the need to iterate step by step over time.
Resource Dynamics
The resource spectrum \(N_R(w)\) is also updated using a similar Euler method, but often with a different equation. For example, if semi-chemostat dynamics are used: \[ \frac{\partial N_R}{\partial t} = r_R(c_R - N_R) - \mu_R N_R. \] The discretised update is: \[ \frac{N_{R,j}^{t+1} - N_{R,j}^t}{\Delta t} = r_R(w_j) (c_R(w_j) - N_{R,j}^{t+1}) - \mu_R(w_j) N_{R,j}^{t+1}. \] Solving for \(N_{R,j}^{t+1}\): \[ N_{R,j}^{t+1} = \frac{N_{R,j}^t + \Delta t r_R(w_j) c_R(w_j)}{1 + \Delta t (r_R(w_j) + \mu_R(w_j))}. \]
Summary of the Algorithm
In each time step \(\Delta t\), the
project() function performs the following steps:
- Calculate Rates: Calculate the biological rates (growth \(g_i\), mortality \(\mu_i\), recruitment \(R_{dd, i}\)) and diffusion coefficients \(d_i\) based on the current population densities \(N^t\).
- Resource Update: Update the resource density \(N_R^{t+1}\) using the explicit rates and implicit density.
-
Consumer Update: Update the consumer densities
\(N_i^{t+1}\) by solving the
tridiagonal system for each species.
- Combine the recruitment flux \(R_{dd, i}\) with the boundary conditions.
- Use a tridiagonal solver (e.g., Thomas algorithm).
- Advance Time: \(t \leftarrow t + \Delta t\).
This split between calculating rates explicitly and solving densities implicitly is what makes the scheme “semi-implicit”.