Lotka-Volterra with harvesting
Lotka-Volterra predation where both species are harvested with a constant effort.
- Name:
lotka_volterra_with_harvesting - Field:
biology/ecology - Formalism: ordinary differential equations
- References: Volterra 1926; Trzaska 2023
Factory
@register_model(
field='biology/ecology',
tags=('predator-prey', 'ode', 'harvesting', 'oscillation'),
references=('Volterra 1926', 'Trzaska 2023'),
defaults={'until': 50.0, 'dt': 0.1},
)
def lotka_volterra_with_harvesting(
prey: float = 16.0,
predators: float = 8.0,
*,
r: float = 1.0,
a: float = 0.1,
b: float = 0.05,
d: float = 0.5,
prey_harvest: float = 0.4,
predator_harvest: float = 0.4,
) -> System:
"""Return Lotka-Volterra predation where both species are harvested with a constant effort.
The populations cycle around `prey = (d + predator_harvest) / b` and `predators = (r - prey_harvest) / a`, which are
also their averages over a cycle. Harvesting raises the average prey and lowers the average predators, and a prey
harvest above `r` drives the predators extinct.
Args:
prey:
The initial number of prey.
predators:
The initial number of predators.
r:
The growth rate of the prey.
a:
The rate at which predation removes prey.
b:
The rate at which predation feeds predators.
d:
The mortality rate of the predators.
prey_harvest:
The harvesting effort on the prey.
predator_harvest:
The harvesting effort on the predators.
Returns:
system:
The model, with the populations `prey` and `predators`.
Examples:
>>> from skcomplex.models.biology.ecology import lotka_volterra_with_harvesting
>>> [interaction.label for interaction in lotka_volterra_with_harvesting().interactions][-2:]
['prey_harvest', 'predator_harvest']
"""
return System(
elements=[PopulationElement('prey', n=prey), PopulationElement('predators', n=predators)],
interactions=[
ExponentialGrowthInteraction('growth', target='prey', r=r),
PredationInteraction('predation', predator='predators', prey='prey', a=a, b=b),
MortalityInteraction('death', target='predators', d=d),
MortalityInteraction('prey_harvest', target='prey', d=prey_harvest),
MortalityInteraction('predator_harvest', target='predators', d=predator_harvest),
],
)
Components
class PopulationElement(CompartmentElement):
"""Population of one species, a compartment whose state `n` holds its size.
Args:
label:
The label of the population.
visual:
How the element is drawn, or `None` for `default_visual_`.
n:
The initial size.
Examples:
>>> from skcomplex.models.biology.common.populations import PopulationElement
>>> PopulationElement('rabbits', n=40).n
40
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual, n=n)
class ExponentialGrowthInteraction(GrowthInteraction):
"""Growth at a constant per capita rate, `dn/dt = r n`, the exponential form of a growth without harvest.
Args:
label:
The label of the interaction.
target:
The path of the element or group whose population grows.
r:
The per capita growth rate, negative for decay.
state:
The name of the state of the target that holds the population.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.common.populations import ExponentialGrowthInteraction
>>> system = System(
... elements=[CompartmentElement('rabbits', n=1.0)],
... interactions=[ExponentialGrowthInteraction('growth', target='rabbits', r=0.5)],
... )
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('rabbits.n')[-1]), 3)
2.718
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
r: Annotated[float, ParameterSpec()] = 1.0,
state: str = 'n',
) -> None:
super().__init__(label, target=target, r=r, form='exponential', state=state)
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the growth `r n`, with no harvest.
Args:
namespace:
The namespace of the symbols the rate can use, with the namespace of the target.
Returns:
rate:
The rate of the population.
"""
return self._get_growth(namespace)
class PredationInteraction(Interaction):
"""Lotka-Volterra predation, where encounters remove prey and feed predators.
The prey loses `a * predator * prey` per unit time and the predator gains `b * predator * prey`.
Args:
label:
The label of the interaction.
predator:
The path of the predator population.
prey:
The path of the prey population.
a:
The attack rate.
b:
The conversion rate of eaten prey into predators.
Examples:
>>> from skcomplex.models.biology.ecology import PredationInteraction
>>> PredationInteraction('hunt', predator='foxes', prey='rabbits', a=0.4, b=0.1).prey
'rabbits'
"""
def __init__(
self,
label: str,
*,
predator: Annotated[str, RoleSpec()],
prey: Annotated[str, RoleSpec()],
a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
) -> None:
super().__init__(label)
self.predator = predator
self.prey = prey
self.a = a
self.b = b
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the rates `-a * predator * prey` and `b * predator * prey`.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
encounters = s.predator.n * s.prey.n
return [RateLaw(s.prey.n, -s.a * encounters), RateLaw(s.predator.n, s.b * encounters)]
class MortalityInteraction(RateInteraction):
"""Death or loss at a constant per capita rate, `dn/dt = -d n`.
Args:
label:
The label of the interaction.
target:
The path of the element or group whose population dies.
d:
The per capita death rate.
state:
The name of the state of the target that holds the population.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.common.populations import MortalityInteraction
>>> system = System(
... elements=[CompartmentElement('foxes', n=8.0)],
... interactions=[MortalityInteraction('death', target='foxes', d=0.5)],
... )
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('foxes.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
state: str = 'n',
) -> None:
super().__init__(label, target=target, state=state)
self.d = d
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the loss `-d n`.
Args:
namespace:
The namespace of the symbols the rate can use, with the namespace of the target.
Returns:
rate:
The rate of the population.
"""
loss: Expression = -namespace.d * getattr(namespace.target, self.state)
return loss
Equations
prey.n: diff(n_prey(t), t) == -a_predation*n_predators(t)*n_prey(t) - d_prey_harvest*n_prey(t) + r_growth*n_prey(t)
predators.n: diff(n_predators(t), t) == b_predation*n_predators(t)*n_prey(t) - d_death*n_predators(t) - d_predator_harvest*n_predators(t)