Case Study 23.1: Porting a Predator-Prey Model from Python to Fortran

"The balance of nature is not a status quo; it is fluid, ever shifting, in a constant state of adjustment." — Rachel Carson, Silent Spring

Executive Summary

A colleague has a working Python script that simulates a predator-prey ecosystem with the classic Lotka-Volterra equations. It is correct but slow: the RHS is called hundreds of thousands of times inside a Python loop, and for the parameter sweep they want to run — thousands of ecosystems, each integrated for centuries of simulated time — pure Python will take hours. Your job is to read their model, port the numerical core to Fortran so it can later be wrapped for Python, and verify the port produces identical numbers on the checkpoints that matter.

This is the everyday reality of scientific computing: you rarely write a solver from a blank page; you inherit a model, understand it, and reimplement its hot core in a faster language. The skill is doing so without changing the science — verifying the port against the original at every step.

Skills applied: - Reading an RHS and reducing a model to the first-order system form (§23.1, §23.4). - Writing a vector RHS as a Fortran rhs_sys function with an assumed-shape argument (§23.4, Chapter 6). - Integrating a system with rk4_sys — one integrator, any RHS (§23.4). - Population-dynamics reasoning: equilibria and conserved quantities as sanity checks (§23.6). - Hand-verifying a port against known checkpoints (the book's no-run discipline).

Background

The Lotka-Volterra equations model two populations — prey $x$ and predator $y$ — that interact: $$x' = a x - b x y, \qquad y' = d x y - c y.$$ Prey grow exponentially ($a x$) but are eaten at a rate proportional to encounters ($b x y$); predators die off ($-c y$) but grow by eating prey ($d x y$). The solutions are famous closed loops: prey rise, predators follow, prey crash, predators starve, prey recover — forever. The system has an equilibrium at $(x^*, y^*) = (c/d,\ a/b)$ where both rates vanish, and it conserves a quantity, so a correct integration keeps the trajectory on a closed curve.

Here is the colleague's Python core, with parameters $a = 1$, $b = 0.1$, $c = 1.5$, $d = 0.075$:

import numpy as np
a, b, c, d = 1.0, 0.1, 1.5, 0.075

def f(t, s):                       # s = [prey, predator]
    x, y = s
    return np.array([a*x - b*x*y, d*x*y - c*y])

def rk4_step(t, s, h):
    k1 = f(t, s)
    k2 = f(t + h/2, s + h/2*k1)
    k3 = f(t + h/2, s + h/2*k2)
    k4 = f(t + h, s + h*k3)
    return s + (h/6)*(k1 + 2*k2 + 2*k3 + k4)

It is a hand-written RK4 on a two-vector — exactly the structure of §23.4. The port is almost mechanical, which is the point: because both languages express the same math, a faithful port is possible and checkable.

Phase 1 — Read the Model and Find Its Fixed Points

Before porting anything, understand the model well enough to test it. The equilibrium is the strongest free check: at $(x^*, y^*) = (c/d, a/b) = (1.5/0.075,\ 1.0/0.1) = (20, 10)$ both derivatives must be exactly zero. Verify by hand: $f_x(20, 10) = 1\cdot20 - 0.1\cdot20\cdot10 = 20 - 20 = 0$, and $f_y(20, 10) = 0.075\cdot20\cdot10 - 1.5\cdot10 = 15 - 15 = 0$. Good — if your port does not return $(0, 0)$ there, it is wrong.

A second, cheaper check: at the off-equilibrium start $(10, 10)$, the prey rate is $1\cdot10 - 0.1\cdot10\cdot10 = 10 - 10 = 0$ (prey exactly balanced for the moment), and the predator rate is $0.075\cdot10\cdot10 - 1.5\cdot10 = 7.5 - 15 = -7.5$ (too few prey, so predators decline). These two numbers — $(0, -7.5)$ — are our port's first target.

Phase 2 — Port the RHS to Fortran

The RHS becomes a Fortran function matching the rhs_sys interface: a scalar time, an assumed-shape state, an array result. The assumed-shape s(:) is what lets the same rk4_sys integrate this two-vector or a million-point PDE unchanged.

module lotka
  implicit none
  integer,  parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: a = 1.0_dp, b = 0.1_dp, c = 1.5_dp, d = 0.075_dp
contains
  function lv(t, s) result(ds)          ! s = [prey, predator]
    real(dp), intent(in) :: t
    real(dp), intent(in) :: s(:)
    real(dp)             :: ds(size(s))
    ds(1) = a*s(1) - b*s(1)*s(2)        ! prey
    ds(2) = d*s(1)*s(2) - c*s(2)        ! predator
  end function lv
end module lotka

Phase 3 — Verify the Checkpoints by Hand

We test the port on the two equilibria rates and two explicit Euler steps — all values we computed by hand in Phase 1, so we can assert them exactly. (Euler, not RK4, for the verification steps, because Euler's arithmetic is short enough to certify on paper.)

program predator_prey
  use lotka
  implicit none
  real(dp) :: s(2), s1(2), s2(2)

  print '(a, 2f11.5)', 'rate at equilibrium (20,10) = ', lv(0.0_dp, [20.0_dp, 10.0_dp])
  print '(a, 2f11.5)', 'rate at start       (10,10) = ', lv(0.0_dp, [10.0_dp, 10.0_dp])

  s  = [10.0_dp, 10.0_dp]
  s1 = s  + 0.1_dp * lv(0.0_dp, s)      ! one explicit Euler step, h = 0.1
  s2 = s1 + 0.1_dp * lv(0.1_dp, s1)     ! a second Euler step
  print '(a, 2f11.5)', 'Euler step 1 (prey,pred)    = ', s1
  print '(a, 2f11.5)', 'Euler step 2 (prey,pred)    = ', s2
end program predator_prey

$ gfortran -std=f2018 -Wall lotka.f90 predator_prey.f90 -o pp && ./pp
rate at equilibrium (20,10) =     0.00000    0.00000
rate at start       (10,10) =     0.00000   -7.50000
Euler step 1 (prey,pred)    =    10.00000    9.25000
Euler step 2 (prey,pred)    =    10.07500    8.55625

Every number is hand-checkable. Step 1: prey stays $10$ (rate $0$), predator falls to $10 + 0.1(-7.5) = 9.25$. Step 2: with prey $10$ and predator $9.25$, the prey rate is $10 - 0.1\cdot10\cdot9.25 = 0.75$ so prey rise to $10.075$; the predator rate is $0.075\cdot10\cdot9.25 - 1.5\cdot9.25 = 6.9375 - 13.875 = -6.9375$ so predator falls to $9.25 - 0.69375 = 8.55625$. The port matches the model exactly. Only now do we trust it with RK4.

Phase 4 — Integrate with RK4 and Watch the Cycle

Swapping the verified Euler diagnostic for rk4_sys(lv, t, s, dt) in a time loop (using the routine from code/example-03-system.f90) gives the production integrator. With $dt = 0.05$ over a few hundred steps, the populations trace the classic cycle. A representative slice of the trajectory (an illustrative Tier-3 table — the exact figures come from the machine, not from hand computation):

The numbers cycle and return — the signature of a conserved system. The diagnostic that proves the integration is faithful is not any single value but the closure of the loop: plot $y$ against $x$ and the trajectory should return to its starting point, not spiral. If it spirals outward, your integrator is injecting energy (plain Euler does exactly this — §23.6); if it spirals inward, it is leaking energy.

There is a sharper, quantitative version of that check. Lotka-Volterra has an exact conserved quantity, $$H(x, y) = d\,x - c\ln x + b\,y - a\ln y,$$ which is constant along every true trajectory. Compute $H$ each step and print its drift: a faithful integrator holds $H$ to within its truncation error, and the drift is your honest error meter even though you have no closed-form solution to compare against. This is the professional move — when you cannot check against the exact answer, check against an invariant the exact answer must obey. RK4 holds $H$ enormously better than Euler, and the drift shrinks by about $16\times$ each time you halve $\Delta t$, the fourth-order signature showing up in a conservation law rather than in a known solution. (You implement this in Extension B.)

Phase 5 — The Payoff and the Handoff

The port is now ready for its real purpose. Wrapped with f2py (Chapter 15), this Fortran lv and its RK4 loop can be called from the colleague's Python analysis code, running the thousands-of-ecosystems sweep at compiled speed while Python still drives the parameter grid and draws the plots. This is the "better together" pattern: Fortran for the hot inner loop, Python for orchestration. The verification we just did — equilibria and two Euler steps — is what lets the colleague trust the faster version: identical checkpoints mean identical science.

Why bother, when SciPy's solve_ivp already integrates Lotka-Volterra in one call? Because the colleague's bottleneck is the number of integrations, not any single one. A parameter sweep of thousands of ecosystems, each marched for hundreds of thousands of steps, calls the RHS billions of times, and every one of those calls is a round trip into the Python interpreter when f is a Python function. Moving f — and the RK4 loop around it — into compiled Fortran removes the interpreter from the inner loop entirely; the Python that remains only launches each run and collects the result. This is the same order-of-magnitude argument the book makes for f2py throughout (Chapter 15): the win is not a faster algorithm but a faster per-call cost, multiplied by billions of calls. The illustrative expectation is a one-to-two-order-of-magnitude speedup for the sweep, with the science bit-for-bit unchanged because the numerical method did not.

Discussion Questions

1. Why is the equilibrium point such a powerful test of a ported RHS? What class of coding errors would it not catch?
2. The verification steps used Euler, but the production run uses RK4. Why verify with the simpler method and only then trust the more accurate one?
3. The trajectory table is labeled "illustrative (Tier 3)." Why can we not honestly present those numbers as hand-computed "expected output," when the Euler-step numbers are presented that way?
4. If plain Euler makes the predator-prey loop spiral outward, what does that tell you physically about what Euler does to a conserved system?

Your Turn: Extensions

• Option A (⭐). Add the equilibrium as a third diagnostic start point: integrate from exactly $(20, 10)$ and confirm the populations stay put (to round-off) for many steps. What step size, if any, makes even the equilibrium drift?
• Option B (⭐⭐). Compute the conserved quantity $H(x, y) = d\,x - c\ln x + b\,y - a\ln y$ at each step and print its drift. Compare Euler versus RK4: how much better does RK4 hold $H$ constant, and how does the drift scale when you halve $dt$?
• Option C (⭐⭐⭐). Add a logistic cap on the prey ($a x(1 - x/K)$ instead of $a x$) and explore how the closed cycles change — do they become a spiral into a stable equilibrium? This is the more realistic model; integrate it and describe the new long-term behavior.

Key Takeaways

• Porting a model means reimplementing the same mathematics in a faster language and proving it with checkpoints — equilibria, conserved quantities, and short hand-traceable steps.
• A vector RHS with an assumed-shape argument plus a generic rk4_sys is all a system of ODEs needs; the predator-prey model differs from the heat solver only in its RHS function.
• For a conserved system, the trajectory's closure, not any single value, is the correctness test — and it exposes whether an integrator adds or leaks energy.
• Verify with the method you can certify by hand (Euler), then trust the method you actually ship (RK4).
• This Fortran core is exactly what Chapter 15 wraps for Python: the fast engine under a Python driver.