Bird swarm with pair potentials
A swarm of self-propelled birds held together by a pair potential that attracts far and repels near.
- Name:
bird_swarm - Field:
biology/collective_behaviour - Formalism: ordinary differential equations
- References: D'Orsogna et al. 2006; Gros 2015
Factory
@register_model(
field='biology/collective_behaviour',
tags=('swarms', 'particles', 'active matter', 'self-propulsion'),
references=("D'Orsogna et al. 2006", 'Gros 2015'),
defaults={'until': 50.0, 'dt': 0.5},
)
def bird_swarm(
n_birds: int = 30,
*,
alpha: float = 1.0,
beta: float = 0.5,
C_a: float = 1.0,
l_a: float = 2.0,
C_r: float = 2.0,
l_r: float = 0.5,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a swarm of self-propelled birds held together by a pair potential that attracts far and repels near.
Every bird is driven to the speed `sqrt(alpha / beta)` and feels the Morse-type potential of every other bird, so
the swarm stays together and organizes into a flock or a rotating mill.
Args:
n_birds:
The number of birds.
alpha:
The propulsion.
beta:
The friction at high speed.
C_a:
The strength of the attraction.
l_a:
The range of the attraction.
C_r:
The strength of the repulsion.
l_r:
The range of the repulsion.
random_state:
The seed or generator of the initial positions and velocities.
Returns:
system:
The model, with the positions `birds.x` and the velocities `birds.x_d1`.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import bird_swarm
>>> bird_swarm(5).elements[0].check().prepare().size_
5
"""
rng = np.random.default_rng(random_state)
positions, velocities = rng.uniform(-3.0, 3.0, (n_birds, 2)), rng.normal(0.0, 1.0, (n_birds, 2))
birds = GroupElement('birds', element=PointElement, size=n_birds, values={'x': positions, 'x_d1': velocities})
potential = MorsePotentialInteraction(
'flocking',
first='birds',
second='birds',
cutoff=10 * l_a,
C_a=C_a,
l_a=l_a,
C_r=C_r,
l_r=l_r,
)
flight = SelfPropulsionInteraction('flight', target='birds', alpha=alpha, beta=beta)
return System(elements=[birds], interactions=[potential, flight], space=ContinuousSpace(n_dim=2))
Components
class PointElement(Element):
"""Point of a continuous space, a member with a vector state `x` of order 2 that locates it.
The laws that act on the point give the second derivative of `x`, and the companion state `x_d1` holds its first
derivative. The boundary of the space acts on `x`.
Args:
label:
The label of the point.
x:
The initial value of the vector state. A number applies to every component.
x_d1:
The initial value of the first derivative of the vector state. A number applies to every component.
visual:
How the point is drawn, or `None` for `default_visual_`.
Attributes:
default_visual_:
Points at the values of `x`.
Examples:
>>> from skcomplex.core.elements import PointElement
>>> from skcomplex.core.spaces import ContinuousSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> system = System(elements=[PointElement('ball', x_d1=[1.0, 0.0])], space=ContinuousSpace(n_dim=2))
>>> Simulation(system, dt=1.0).check().prepare().run(until=2.0).trajectory_.values('ball.x')[-1].tolist()
[2.0, 0.0]
"""
def __init__(
self,
label: str,
x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2, role='position')] = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.x_d1 = x_d1
def _prepare(self, context: BaseSpace | None = None) -> None:
"""Set `default_visual_`, points at the values of `x`.
Args:
context:
The space the element inherits from its system, which the element does not read.
"""
self.default_visual_ = VisualSpec('points', {'position': 'x'})
class MorsePotentialInteraction(PairCouplingInteraction):
"""Attraction at long range and repulsion at short range, `C_r exp(-r / l_r) - C_a exp(-r / l_a)`.
Args:
label:
The label of the interaction.
first:
The path of the group or particle of the first members.
second:
The path of the group or particle of the second members.
cutoff:
The distance beyond which pairs do not interact.
position:
The name of the position state of the members.
period:
The side of a periodic box, for minimum image distances, or `None` for an open space.
C_a:
The strength of the attraction.
l_a:
The range of the attraction.
C_r:
The strength of the repulsion.
l_r:
The range of the repulsion.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import MorsePotentialInteraction
>>> MorsePotentialInteraction('flocking', first='birds', second='birds', cutoff=10.0).l_a
2.0
"""
def __init__(
self,
label: str,
first: str,
second: str,
cutoff: float = 1.0,
position: str = 'x',
*,
period: float | None = None,
C_a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
l_a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
C_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
l_r: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
) -> None:
super().__init__(
label,
first,
second,
cutoff=cutoff,
position=position,
potential=self._get_pair_energy,
period=period,
)
self.C_a = C_a
self.l_a = l_a
self.C_r = C_r
self.l_r = l_r
@staticmethod
def _get_pair_energy(s: SimpleNamespace, r: Expression) -> Expression:
"""Return `C_r exp(-r / l_r) - C_a exp(-r / l_a)`."""
return s.C_r * exp(-r / s.l_r) - s.C_a * exp(-r / s.l_a)
class SelfPropulsionInteraction(RateInteraction):
"""Rayleigh self-propulsion, `v' = (alpha - beta |v|²) v`, which drives points to the speed `sqrt(alpha / beta)`.
It is a rate on the position `x`, a state of order 2, whose velocity `v` is the first derivative `x_d1`.
Args:
label:
The label of the interaction.
target:
The path of the particle or group.
alpha:
The propulsion.
beta:
The friction at high speed.
Examples:
>>> from skcomplex.models.biology.collective_behaviour import SelfPropulsionInteraction
>>> SelfPropulsionInteraction('flight', target='birds', alpha=1.5).alpha
1.5
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
) -> None:
super().__init__(label, target=target, state='x')
self.alpha = alpha
self.beta = beta
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the velocity change `(alpha - beta |v|²) v`.
Args:
namespace:
The namespace of the symbols the rate can use, with the namespace of the target.
Returns:
rate:
The change of the velocity.
"""
velocity = namespace.target.x_d1
return (namespace.alpha - namespace.beta * velocity.dot_product(velocity)) * velocity