Predator-Prey Population Dynamics: Symbolic ODEs with Interval Uncertainty

Overview & Architecture

TheThe Lotka-Volterra model describes the two populations with a pair of coupled ordinary differential equations: \( \frac{dR}{dt} = aR - bRF \) (rabbits grow at rate a, get eaten at a rate proportional to how often a rabbit and a fox meet) and \( \frac{dF}{dt} = -cF + dRF \) (foxes starve at rate c, breed at a rate proportional to that same meeting frequency). self_hosting/lib/symbolic_ode.patlang integrates this two ways, deliberately kept separate rather than blended into one method:

• A plain fourth-order Runge-Kutta step (lv_step_rk4) — ordinary point-value numerical integration, no intervals at all. This is the accurate, practical method, and it produces the population-over-time and fox-versus-rabbit charts below.
• An interval Euler step (lv_step_interval), built entirely from the sign-flip-aware sym_interval_mul/sym_interval_add already proven on the Symbolic Math page. Given a starting population that's only known to within a small range rather than exactly, this propagates that uncertainty through the real nonlinear coupling term R·F at every step.

The second one needs a real caveat, stated plainly rather than left for a reader to discover by disappointment: this is uncertainty propagation, not a rigorous enclosure the way the sqrt(2) bracket on the Interval Arithmetic page is. A genuinely rigorous interval ODE solver needs an a priori bound on the solution across each step to certify the local truncation error — the interval Taylor and interval Euler methods Moore's own book covers¹ — which this does not attempt. What it does show for real is the wrapping effect: watch the shaded band below widen far faster than the actual uncertainty would, until it crosses into physically meaningless negative populations. That failure mode is exactly why rigorous interval ODE solving is a substantially harder problem than rigorous interval arithmetic on its own, not a bug in this demo.

Run It

Run simulation

(not run yet)

(not run yet)

Parameters

FixedFixed for this demo: a = 1.0, b = 0.1, c = 1.5, d = 0.075, starting at 10 rabbits and 5 foxes — ordinary illustrative values chosen to produce a clean, clearly-oscillating cycle (period roughly 5 time units), not fitted to any real population data. The uncertainty-propagation charts start from a ±1% range on each population (R₀ ∈ [9.9, 10.1], F₀ ∈ [4.9, 5.1]) and are deliberately cut off at t = 1.8, once the band has already crossed into negative populations — continuing further only makes the chart's own scale useless, not the demonstration more honest.

# Lotka-Volterra predator-prey ODE integration
# (self_hosting/lib/symbolic_ode.patlang, built on
# self_hosting/lib/symbolic.patlang's Interval type).
#
# dR/dt = a*R - b*R*F
# dF/dt = -c*F + d*R*F
#
# Two integrators, deliberately not one: lv_step_rk4 is a plain point-value
# fourth-order Runge-Kutta step -- the accurate practical method. lv_step_
# interval is an interval Euler step built from sym_interval_add/mul --
# uncertainty PROPAGATION through the real R*F coupling term, not a
# rigorous ENCLOSURE (see the page this drives for why those differ).

make a function called lv_deriv takes r, f, a, b, c, d returns pair
  let dr = (a * r) - (b * r * f)
  let df = ((0 - c) * f) + (d * r * f)
  return [dr, df]
end

make a function called lv_step_rk4 takes r, f, a, b, c, d, h returns pair
  let k1 = lv_deriv(r, f, a, b, c, d)
  let k2 = lv_deriv(r + ((h / 2) * k1[0]), f + ((h / 2) * k1[1]), a, b, c, d)
  let k3 = lv_deriv(r + ((h / 2) * k2[0]), f + ((h / 2) * k2[1]), a, b, c, d)
  let k4 = lv_deriv(r + (h * k3[0]), f + (h * k3[1]), a, b, c, d)
  let r2 = r + ((h / 6) * (k1[0] + (2 * k2[0]) + (2 * k3[0]) + k4[0]))
  let f2 = f + ((h / 6) * (k1[1] + (2 * k2[1]) + (2 * k3[1]) + k4[1]))
  return [r2, f2]
end

make a function called lv_deriv_interval takes r, f, a, b, c, d returns pair
  let rf = sym_interval_mul(r, f)
  let dr = sym_interval_sub(sym_interval_mul(sym_interval(a, a), r), sym_interval_mul(sym_interval(b, b), rf))
  let df = sym_interval_add(sym_interval_mul(sym_interval(0 - c, 0 - c), f), sym_interval_mul(sym_interval(d, d), rf))
  return [dr, df]
end

make a function called lv_step_interval takes r, f, a, b, c, d, h returns pair
  let d1 = lv_deriv_interval(r, f, a, b, c, d)
  let r2 = sym_interval_add(r, sym_interval_mul(sym_interval(h, h), d1[0]))
  let f2 = sym_interval_add(f, sym_interval_mul(sym_interval(h, h), d1[1]))
  return [r2, f2]
end

let a = 1.0
let b = 0.1
let c = 1.5
let d = 0.075

print("RK4")
let h1 = 0.05
let n1 = 500
let r = 10.0
let f = 5.0
let i = 0
while i <= n1 do
  let t = i * h1
  print(t + "," + r + "," + f)
  let step = lv_step_rk4(r, f, a, b, c, d, h1)
  let r = step[0]
  let f = step[1]
  let i = i + 1
end

print("INTERVAL")
let h2 = 0.01
let n2 = 180
let ri = sym_interval(9.9, 10.1)
let fi = sym_interval(4.9, 5.1)
let i = 0
while i <= n2 do
  let t = i * h2
  print(t + "," + sym_ival_lo(ri) + "," + sym_ival_hi(ri) + "," + sym_ival_lo(fi) + "," + sym_ival_hi(fi))
  let stepi = lv_step_interval(ri, fi, a, b, c, d, h2)
  let ri = stepi[0]
  let fi = stepi[1]
  let i = i + 1
end
print("")

References