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.
Summary of simulation algorithm
In each time step \(\Delta t\), the project() function performs the following steps:
-
Calculate Rates: Calculate the biological rates (growth \(g\), mortality \(\mu\), recruitment \(R_{dd}\)) and diffusion coefficients \(d\) 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^{t+1}\) by solving the tridiagonal system for each species.
- Combine the recruitment flux \(R_{dd}\) 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”.
The weight grid
The model dynamics are described by partial differential equations (PDEs) for the consumer number density \(N(w)\) of each species 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. In this vignette we suppress the species index and write the equations for one species; the same scheme is applied to each species.
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}. \tag{1}\] 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). \tag{2}\] Note that \(\Delta w_j\) increases with \(w_j\), which is appropriate for a logarithmic grid.
The number density \(N(w)\) is approximated by a constant value \(N_j\) within each bin. Thus \(N_j\) represents the density at the start of the bin \(w_j\). The number of individuals in the \(j\)-th bin is \(N_j \Delta w_j\).
The Transport Equation
The time evolution of the number density \(N(w)\) is described by the McKendrick-von Foerster equation with an added diffusion term: \[
\frac{\partial N}{\partial t} + \frac{\partial}{\partial w} \left( g N - \frac{1}{2}\frac{\partial(d N)}{\partial w} \right) = -\mu N
\tag{3}\] where \(g(w)\) is the somatic growth rate, \(d(w)\) is the diffusion coefficient and \(\mu(w)\) is the total mortality rate.
We discretise this equation using a finite volume scheme. On the grid we write \(g_j = g(w_j)\), \(d_j = d(w_j)\), and \(\mu_j = \mu(w_j)\).
Discretisation of the Fluxes
The term inside the derivative with respect to \(w\) is the flux \(J(w)\): \[ J(w) = g(w) N(w) - \frac{1}{2}\frac{\partial(d(w) N(w))}{\partial w} \tag{4}\] We consider the \(j\)-th size bin \([w_j, w_{j+1}]\). Integrating the conservation equation over this bin gives: \[ \int_{w_j}^{w_{j+1}} \frac{\partial N_j}{\partial t}dw + J(w_{j+1}) - J(w_j) = -\int_{w_j}^{w_{j+1}} \mu(w) N(w) dw \tag{5}\] Approximating the integral and dividing by \(\Delta w_j\): \[ \frac{\partial N_j}{\partial t}+ \frac{J_{j+1} - J_j}{\Delta w_j} = -\mu_j N_j \tag{6}\] where \(J_j\) represents the flux at the boundary \(w_j\).
Advective Flux
For the advective part of the flux \(g N\), we use an upwind approximation. Since fish grow from smaller to larger sizes (\(g > 0\)), the flux at boundary \(w_j\) is determined by the density in the bin below (\(j-1\)): \[ J_j^{adv} = g_{j-1} N_{j-1} \tag{7}\] (Note: for \(j=j_{min}\), this boundary flux is the recruitment \(R_{dd}\)).
Diffusive Flux
For the diffusive part of the flux \(-\frac{1}{2}\frac{\partial(d N)}{\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_j^{diff} \approx -\frac{1}{2} \frac{(d N)_j - (d N)_{j-1}}{w_j - w_{j-1}} = -\frac{1}{2} \frac{d_j N_j - d_{j-1} N_{j-1}}{\Delta w_{j-1}} \tag{8}\] Note the use of \(N_j\) (density in bin \(j\)) and \(N_{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_j = g_{j-1} N_{j-1} - \frac{1}{2} \frac{d_j N_j - d_{j-1} N_{j-1}}{\Delta w_{j-1}} \tag{9}\] And similarly at boundary \(w_{j+1}\): \[ J_{j+1} = g_j N_j - \frac{1}{2} \frac{d_{j+1} N_{j+1} - d_j N_j}{\Delta w_j} \tag{10}\]
Semi-Implicit Time Discretisation
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\) are evaluated at time \(t+1\), but the rates (\(g\), \(\mu\), \(d\)) are evaluated at time \(t\). Using a fully implicit scheme would require solving a nonlinear system at each time step, which is more computationally expensive. Thus, the discretised equation becomes: \[
\frac{N_j^{t+1} - N_j^t}{\Delta t} + \frac{1}{\Delta w_j} \left( J_{j+1}^{t+1} - J_j^{t+1} \right) = -\mu_j N_j^{t+1}
\tag{11}\] where the fluxes \(J_j^{t+1}\) are calculated using the densities at time \(t+1\) but the rates at time \(t\). To simplify the notation we will drop the explicit time indices on the rates, but it is important to remember that they are evaluated at time \(t\). The fluxes at time \(t+1\) are: \[
\begin{aligned}
J_{j+1}^{t+1} &= g_j N_j^{t+1} - \frac{1}{2} \frac{d_{j+1} N_{j+1}^{t+1} - d_j N_j^{t+1}}{\Delta w_j} \\
J_j^{t+1} &= g_{j-1} N_{j-1}^{t+1} - \frac{1}{2} \frac{d_j N_j^{t+1} - d_{j-1} N_{j-1}^{t+1}}{\Delta w_{j-1}}.
\end{aligned}
\tag{12}\] This leads to a linear system of the form: \[
A_j N_{j-1}^{t+1} + B_j N_j^{t+1} + C_j N_{j+1}^{t+1} = S_j
\tag{13}\] where \(S_j = N_j^t\). This is a tridiagonal system for each species, which can be solved efficiently (e.g., using the Thomas algorithm).
The coefficients are: \[
\begin{aligned}
A_j &= -\frac{\Delta t}{\Delta w_j} \left( g_{j-1} + \frac{1}{2} \frac{d_{j-1}}{\Delta w_{j-1}} \right) \\
C_j &= -\frac{\Delta t}{\Delta w_j} \left( \frac{1}{2} \frac{d_{j+1}}{\Delta w_j} \right) \\
B_j &= 1 + \Delta t \mu_j + \frac{\Delta t}{\Delta w_j} \left( g_j + \frac{1}{2} \frac{d_j}{\Delta w_j} + \frac{1}{2} \frac{d_j}{\Delta w_{j-1}} \right)
\end{aligned}
\tag{14}\]
Boundary Conditions
At the smallest size (\(j=j_{min}\)): The flux entering the grid is determined by recruitment: \[ J_{j_{min}}^{t+1} = R_{dd} \tag{15}\] (We assume diffusive flux at the lower boundary is negligible or incorporated into \(R_{dd}\)). The equation for the first bin becomes: \[ \frac{N_{j_{min}}^{t+1} - N_{j_{min}}^t}{\Delta t} + \frac{J_{j_{min}+1}^{t+1} - R_{dd}}{\Delta w_{j_{min}}} = -\mu_{j_{min}} N_{j_{min}}^{t+1} \tag{16}\] This involves \(N_{j_{min}}^{t+1}\) and \(N_{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_{j_{min}} = 0\).
- The upward diffusion term from below the boundary is omitted, so \(B_{j_{min}}\) does not have the \(\frac{1}{2} \frac{d_{j_{min}}}{\Delta w_{j_{min}-1}}\) component: \[ B_{j_{min}} = 1 + \Delta t \mu_{j_{min}} + \frac{\Delta t}{\Delta w_{j_{min}}} \left( g_{j_{min}} + \frac{1}{2} \frac{d_{j_{min}}}{\Delta w_{j_{min}}} \right) \tag{17}\]
- The coefficient \(C_{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_{j_{min}}\): \[ S_{j_{min}} = N_{j_{min}}^t + \frac{\Delta t}{\Delta w_{j_{min}}} R_{dd} \tag{18}\]
(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_{j_{max}+1} = 0\). The flux leaving the grid is: \[ J_{j_{max}+1}^{t+1} = g_{j_{max}} N_{j_{max}}^{t+1} - \frac{1}{2} \frac{0 - d_{j_{max}} N_{j_{max}}^{t+1}}{\Delta w_{j_{max}}} \tag{19}\] This means the term \(C_{j_{max}}\) multiplying \(N_{j_{max}+1}^{t+1}\) is not needed, so \(C_{j_{max}} = 0\). The coefficients \(A_{j_{max}}\) and \(B_{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 \tag{20}\] 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} \tag{21}\] The coefficient of the second derivative represents the numerical diffusivity: \[ D_{num} = \frac{g \Delta w}{2} (1 + C) \tag{22}\] 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)) \tag{23}\] 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 \tag{24}\] 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\)).
Order of Accuracy
The scheme is first order in time. To see this, expand the backward difference around the new time level: \[
\frac{N_j^{t+1} - N_j^t}{\Delta t}
= \frac{\partial N_j}{\partial t}(t+\Delta t)
- \frac{\Delta t}{2}\frac{\partial^2 N_j}{\partial t^2}(t+\Delta t)
+ O(\Delta t^2).
\tag{25}\] Thus the time discretisation has a truncation error of order \(O(\Delta t)\). Evaluating the rates \(g\), \(d\) and \(\mu\) at time \(t\) instead of \(t+1\) also introduces only an \(O(\Delta t)\) error, provided the rates vary smoothly in time. The semi-implicit scheme is therefore first order in \(\Delta t\).
For the spatial discretisation, the upwind advective flux is the limiting term. For smooth \(gN\), \[
\frac{(gN)_j - (gN)_{j-1}}{\Delta w_{j-1}}
= \frac{\partial(gN)}{\partial w}(w_j)
- \frac{\Delta w_{j-1}}{2}\frac{\partial^2(gN)}{\partial w^2}(w_j)
+ O(\Delta w_{j-1}^2),
\tag{26}\] so the advective part has an \(O(\Delta w_j)\) spatial truncation error. The diffusive part uses a centred difference for the flux and is of higher order on a locally uniform grid, but the overall spatial order is still set by the upwind advective term.
On the logarithmic grid, writing \(\beta = w_{j+1}/w_j\), we have \[
\Delta w_j = w_j(\beta - 1) = w_j\left(\log\beta + O((\log\beta)^2)\right).
\tag{27}\] Thus, at fixed body size \(w_j\), first order in \(\Delta w_j\) is equivalently first order in \(\log\beta\). Combining the time and space errors, the scheme is \[
O(\Delta t) + O(\Delta w_j)
\tag{28}\] locally, or \(O(\Delta t + \log\beta)\) on the logarithmic grid.
Predictor-Corrector Time Stepping
A straightforward way to make the time discretisation second order is to combine a predictor step for the rates with a Crank-Nicolson corrector for the densities. The aim is to keep the useful property that, once the rates are fixed, the density update is still a tridiagonal linear solve.
Let \(L(r)\) denote the spatial operator for the consumer spectrum when the rates \[
r = (g, d, \mu)
\tag{29}\] are fixed, and let \(q(R_{dd})\) denote the recruitment source at the lower boundary. The semi-implicit Euler method used above can be written schematically as \[
\frac{N^{t+1} - N^t}{\Delta t} = L(r^t) N^{t+1} + q(R_{dd}^t).
\tag{30}\] This is first order because the rates and recruitment are only evaluated at the start of the time step.
The predictor-corrector method proceeds as follows:
-
Predict: Use the existing semi-implicit Euler scheme with rates \(r^t\) to get a provisional value \(\hat{N}^{t+1}\).
-
Recalculate rates: Evaluate provisional end-of-step rates \(\hat{r}^{t+1}\) and recruitment \(\hat{R}_{dd}^{t+1}\) from \(\hat{N}^{t+1}\).
-
Approximate midpoint rates: Use \[
r^{t+1/2} = \frac{1}{2}\left(r^t + \hat{r}^{t+1}\right),
\tag{31}\] so in particular \[
g_j^{t+1/2} = \frac{1}{2}\left(g_j^t + \hat{g}_j^{t+1}\right), \quad
d_j^{t+1/2} = \frac{1}{2}\left(d_j^t + \hat{d}_j^{t+1}\right), \quad
\mu_j^{t+1/2} = \frac{1}{2}\left(\mu_j^t + \hat{\mu}_j^{t+1}\right).
\tag{32}\] The recruitment flux is averaged in the same way: \[
R_{dd}^{t+1/2} = \frac{1}{2}\left(R_{dd}^t + \hat{R}_{dd}^{t+1}\right).
\tag{33}\]
-
Correct: Do a Crank-Nicolson update using the midpoint rates: \[
\frac{N^{t+1} - N^t}{\Delta t}
= \frac{1}{2}\left(L(r^{t+1/2})N^t + L(r^{t+1/2})N^{t+1}\right)
+ q(R_{dd}^{t+1/2}).
\tag{34}\]
In flux form, the corrector equation for an interior bin is \[
\begin{aligned}
\frac{N_j^{t+1} - N_j^t}{\Delta t}
&+ \frac{1}{2\Delta w_j}
\left[
\left(J_{j+1}^{t+1} - J_j^{t+1}\right)
+ \left(J_{j+1}^t - J_j^t\right)
\right] \\
&= -\frac{1}{2}\mu_j^{t+1/2}\left(N_j^{t+1} + N_j^t\right),
\end{aligned}
\tag{35}\] where all fluxes use the midpoint rates. For example, \[
J_j^{t+1}
= g_{j-1}^{t+1/2} N_{j-1}^{t+1}
- \frac{1}{2}\frac{
d_j^{t+1/2} N_j^{t+1}
- d_{j-1}^{t+1/2} N_{j-1}^{t+1}
}{\Delta w_{j-1}},
\tag{36}\] and \(J_j^t\) is the same expression with \(N^t\) instead of \(N^{t+1}\). The lower boundary uses \[
J_{j_{min}}^{t+1} = J_{j_{min}}^t = R_{dd}^{t+1/2}.
\tag{37}\]
Because the midpoint rates are fixed during the corrector step, the unknowns still enter only through \(N_{j-1}^{t+1}\), \(N_j^{t+1}\) and \(N_{j+1}^{t+1}\). The corrector is therefore again a tridiagonal system. The difference from the first-order scheme is that the old-time fluxes and mortality terms are moved to the right-hand side, while the new-time terms supply the tridiagonal matrix.
If the provisional predictor has the usual one-step accuracy, the averaged rates approximate the true midpoint rates to second order. The Crank-Nicolson corrector is then second order in \(\Delta t\) for smooth solutions: \[
\text{time error} = O(\Delta t^2).
\tag{38}\] This does not change the spatial order of the scheme, which remains first order because of the upwind advective flux. The combined order would therefore be \[
O(\Delta t^2) + O(\Delta w_j),
\tag{39}\] or \(O(\Delta t^2 + \log\beta)\) on the logarithmic grid.
There are two practical cautions. First, this doubles the number of rate evaluations per time step. Second, Crank-Nicolson-type methods are less damping than backward Euler, so for large time steps they can show oscillations even when they are formally stable. Positivity and robustness would therefore need to be checked carefully before using this as the default time stepper.
Accuracy comparison
This vignette compares the first-order consumer density update, selected with method = "euler", with the predictor-corrector update selected with method = "predictor_corrector", on the North Sea example model NS_params.
The comparison is intentionally modest. It is meant as a reproducible smoke test for speed and time-step sensitivity, not as a comprehensive benchmark.
We compare final consumer spectra against a smaller-time-step reference solution from the predictor-corrector method. Because both methods use the same spatial discretisation, this test is only about the time discretisation error on the fixed NS_params weight grid.
Code
t_max <- 8
dt_values <- c(0.4, 0.2, 0.1, 0.05, 0.025)
reference_dt <- 0.4 / 2^6
params <- NS_params
initial_effort(params) <- 4
relative_l2_error <- function(x, reference) {
x_final <- finalN(x)
reference_final <- finalN(reference)
sqrt(sum((x_final - reference_final)^2)) / sqrt(sum(reference_final^2))
}
reference <- project(
params,
dt = reference_dt,
t_max = t_max,
method = "predictor_corrector"
)
accuracy <- do.call(rbind, lapply(dt_values, function(dt) {
euler <- project(params, dt = dt, t_max = t_max, t_save = t_max,
method = "euler")
pc <- project(params, dt = dt, t_max = t_max, t_save = t_max,
method = "predictor_corrector")
data.frame(
dt = dt,
euler_error = relative_l2_error(euler, reference),
predictor_corrector_error = relative_l2_error(pc, reference)
)
}))
accuracy
dt euler_error predictor_corrector_error
1 0.400 0.06986846 0.2499180214
2 0.200 0.07985741 0.0312481139
3 0.100 0.06201718 0.0096656723
4 0.050 0.03985931 0.0029264864
5 0.025 0.02282481 0.0008933327
Code
plot(
euler_error ~ dt,
data = accuracy,
log = "xy",
type = "b",
pch = 16,
xlab = expression(Delta * t),
ylab = "relative error in final consumer spectrum",
ylim = range(accuracy$euler_error, accuracy$predictor_corrector_error)
)
lines(
predictor_corrector_error ~ dt,
data = accuracy,
type = "b",
pch = 17,
col = 2
)
legend(
"bottomright",
legend = c('method = "euler"', 'method = "predictor_corrector"'),
pch = c(16, 17),
col = c(1, 2),
lty = 1,
bty = "n"
)
Speed
The second-order method does roughly twice the work per time step. In the full predictor-corrector form it performs one predictor solve and then recalculates the rates before doing the corrector solve.
Interpretation
On short runs the predictor-corrector method may not always show a clean second-order convergence slope because the comparison is affected by the fixed spatial grid, the first-order resource update, nonlinear rate feedbacks and ordinary timing noise.
Because the predictor-corrector method is more expensive per time step, the relevant comparison is between the Euler method at a given dt and the predictor-corrector method at twice the dt so that both have similar run times.
Numerical Diffusion for the Predictor-Corrector Method
The predictor-corrector method uses the same upwind approximation in size as the first-order method, so it still has numerical diffusion. However, the Crank-Nicolson corrector changes the time-stepping contribution to that diffusion.
To see this, again consider the advection-only problem with constant growth rate \(g\) on a locally uniform grid with spacing \(\Delta w\): \[
\frac{\partial N}{\partial t} + g\frac{\partial N}{\partial w} = 0.
\tag{40}\] With fixed rates, the corrector step reduces to the Crank-Nicolson upwind scheme \[
\frac{N_j^{t+1} - N_j^t}{\Delta t}
+ \frac{g}{2}
\left[
\frac{N_j^{t+1} - N_{j-1}^{t+1}}{\Delta w}
+ \frac{N_j^t - N_{j-1}^t}{\Delta w}
\right] = 0.
\tag{41}\]
Expand this equation about the midpoint \((w_j, t + \Delta t / 2)\). The centred time difference gives \[
\frac{N_j^{t+1} - N_j^t}{\Delta t}
= \frac{\partial N}{\partial t}
+ O(\Delta t^2).
\tag{42}\] The average of the two upwind spatial differences gives \[
\frac{1}{2}
\left[
\frac{N_j^{t+1} - N_{j-1}^{t+1}}{\Delta w}
+ \frac{N_j^t - N_{j-1}^t}{\Delta w}
\right]
= \frac{\partial N}{\partial w}
- \frac{\Delta w}{2}\frac{\partial^2 N}{\partial w^2}
+ O(\Delta w^2) + O(\Delta t^2).
\tag{43}\] Substituting these expansions into the scheme gives the modified equation \[
\frac{\partial N}{\partial t} + g\frac{\partial N}{\partial w}
= \frac{g\Delta w}{2}\frac{\partial^2 N}{\partial w^2}
+ O(\Delta w^2) + O(\Delta t^2).
\tag{44}\] Thus the numerical diffusivity in the usual advection-diffusion form is \[
D_{num}^{PC} = \frac{g\Delta w}{2}.
\tag{45}\] In the mizer notation, where the diffusion term is written with a factor \(1/2\) inside the flux, this corresponds to the effective diffusion parameter \[
d_{num}^{PC}(w) \approx g(w)\Delta w.
\tag{46}\] On the logarithmic grid this becomes \[
d_{num}^{PC}(w) \approx g(w)w\log\beta.
\tag{47}\]
Compared with the first-order semi-implicit scheme, \[
d_{num}^{Euler}(w) \approx g(w)\Delta w + g(w)^2\Delta t,
\tag{48}\] the predictor-corrector method removes the leading artificial diffusion proportional to \(\Delta t\). The remaining numerical diffusion is the spatial upwind diffusion, proportional to \(\Delta w\) or equivalently to \(\log\beta\) on the logarithmic grid. This is why decreasing \(\Delta t\) eventually stops improving the error much: once the time-stepping error is small, the spatial upwind error and other fixed discretisation errors dominate.
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_j^{t+1} = N_j^t = N_j^*\).
Substituting this into our discretised linear system: \[
A_j N_{j-1}^* + B_j N_j^* + C_j N_{j+1}^* = S_j
\tag{49}\] Recall that for \(j > j_{min}\), \(S_j = N_j^t\). The equation simplifies to: \[
A_j N_{j-1}^* + (B_j - 1) N_j^* + C_j N_{j+1}^* = 0
\tag{50}\]
To find the steady-state population densities \(N^*\), we formulate a new time-independent tridiagonal system: \[
\tilde{A}_j N_{j-1}^* + \tilde{B}_j N_j^* + \tilde{C}_j N_{j+1}^* = \tilde{S}_j
\tag{51}\] 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}_j &= \frac{A_j}{\Delta t} = -\frac{1}{\Delta w_j} \left( g_{j-1} + \frac{1}{2} \frac{d_{j-1}}{\Delta w_{j-1}} \right) \\
\tilde{C}_j &= \frac{C_j}{\Delta t} = -\frac{1}{\Delta w_j} \left( \frac{1}{2} \frac{d_{j+1}}{\Delta w_j} \right) \\
\tilde{B}_j &= \frac{B_j - 1}{\Delta t} = \mu_j + \frac{1}{\Delta w_j} \left( g_j + \frac{1}{2} \frac{d_j}{\Delta w_j} + \frac{1}{2} \frac{d_j}{\Delta w_{j-1}} \right)
\end{aligned}
\tag{52}\] Notice that \(\tilde{A}_j\) and \(\tilde{C}_j\) are exactly the expressions for \(A_j\) and \(C_j\) evaluated at \(\Delta t = 1\). Similarly, \(\tilde{B}_j\) is exactly the expression for \(B_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_{j_{min}} = N_{j_{min}}^t + \frac{\Delta t}{\Delta w_{j_{min}}} R_{dd}
\tag{53}\] Following the same logic of setting \(N^{t+1} = N^t = N^*\) and dividing by \(\Delta t\), the right-hand side vector \(\tilde{S}_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}\) when evaluated at \(\Delta t = 1\), we get our new source vector: \[
\tilde{S}_{j_{min}} = \frac{R_{dd}}{\Delta w_{j_{min}}}
\tag{54}\] For all other \(j > j_{min}\), \(\tilde{S}_j = 0\).
The boundary condition modifications at the edges of the grid remain the same conceptually: \(\tilde{A}_{j_{min}} = 0\), the upward diffusion term is omitted from \(\tilde{B}_{j_{min}}\), and \(\tilde{C}_{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_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.
Steady-State Solution for the Predictor-Corrector Method
The predictor-corrector method changes the time stepping, but it does not change the steady state that is obtained when the rates are evaluated at that steady state. To see this, write the Crank-Nicolson corrector with fixed midpoint rates as \[
\frac{N^{t+1} - N^t}{\Delta t}
= \frac{1}{2}\left(L N^t + L N^{t+1}\right) + q,
\tag{55}\] where \(L\) is the spatial transport-and-mortality operator built from the midpoint rates and \(q\) is the recruitment source at the lower boundary. At a steady state, \(N^{t+1} = N^t = N^*\), so this becomes simply \[
0 = L N^* + q.
\tag{56}\]
Thus the predictor-corrector method has the same fixed-point equation as the first-order method. The predictor step also becomes irrelevant at the fixed point: if \(N^t = N^*\), then the predicted \(\hat{N}^{t+1}\) is also \(N^*\) up to the residual of the steady-state equation, so \[
r^t = \hat{r}^{t+1} = r^{t+1/2}, \qquad
R_{dd}^t = \hat{R}_{dd}^{t+1} = R_{dd}^{t+1/2}.
\tag{57}\] The rates used in the steady-state calculation are therefore just the rates evaluated at \(N^*\). The predictor-corrector method affects the transient path to the steady state, not the steady state itself.
