Cucker-Smale flocking
Cucker-Smale flocking, where every agent aligns its velocity with all others, weighted by distance.
- Name:
cucker_smale_flocking - Field:
physics/active_matter - Formalism: ordinary differential equations
- References: Cucker and Smale 2007; Vicsek and Zafeiris 2012; Bellomo et al. 2017
Factory
@register_model(
field='physics/active_matter',
tags=('flocking', 'ode', 'alignment', 'pair interaction'),
references=('Cucker and Smale 2007', 'Vicsek and Zafeiris 2012', 'Bellomo et al. 2017'),
defaults={'until': 20.0, 'dt': 0.5, 'solver': RK4Solver(dt=0.05)},
)
def cucker_smale_flocking(
n_agents: int = 50,
*,
K: float = 1.0,
beta: float = 0.25,
size: float = 10.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return Cucker-Smale flocking, where every agent aligns its velocity with all others, weighted by distance.
The mean velocity is conserved, and for `beta <= 1/2` the velocities converge to it from any start.
Args:
n_agents:
The number of agents.
K:
The coupling strength.
beta:
The decay exponent of the communication weight.
size:
The side of the square in which the agents start.
random_state:
The seed or generator of the initial positions and velocities.
Returns:
system:
The model, with the positions `birds.x` and velocities `birds.x_d1`.
Examples:
>>> from skcomplex.models.physics.active_matter import cucker_smale_flocking
>>> cucker_smale_flocking(5).interactions[0].n_agents
5
"""
rng = np.random.default_rng(random_state)
birds = GroupElement(
'birds',
element=PointElement,
size=n_agents,
values={'x': rng.uniform(0.0, size, (n_agents, 2)), 'x_d1': rng.normal(0.0, 1.0, (n_agents, 2))},
visual=VisualSpec('points', {'position': 'x', 'color_by': 'x_d1', 'states': ()}),
)
alignment = CuckerSmaleAlignmentInteraction(
'align',
first='birds',
second='birds',
cutoff=1e6,
K=K,
beta=beta,
n_agents=n_agents,
)
return System(elements=[birds], interactions=[alignment], 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 CuckerSmaleAlignmentInteraction(PairCouplingInteraction):
"""Cucker-Smale alignment, `x_i'' += K / N psi(r_ij) (x_j' - x_i')` with `psi(r) = (1 + r^2)^(-beta)`.
Every pair of members within the cutoff pulls their velocities, the first derivatives of their positions of order
2, towards each other with a weight that decays with their distance.
Args:
label:
The label of the interaction, unique among the parts of its system.
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 align, large for the all-to-all model.
position:
The name of the vector state of order 2 that holds the positions of the members.
K:
The strength of the alignment.
beta:
The decay exponent of the weight.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
n_agents:
The number of agents `N` that normalizes the alignment.
Examples:
>>> from skcomplex.models.physics.active_matter import CuckerSmaleAlignmentInteraction
>>> CuckerSmaleAlignmentInteraction('align', first='birds', second='birds', beta=0.25).beta
0.25
"""
def __init__(
self,
label: str,
first: str,
second: str,
cutoff: float = 1.0,
position: str = 'x',
*,
K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.25,
period: float | None = None,
n_agents: int = 1,
) -> None:
super().__init__(
label,
first,
second,
cutoff=cutoff,
position=position,
write=position,
function=self._align,
period=period,
)
self.K = K
self.beta = beta
self.n_agents = n_agents
def _align(self, namespace: SimpleNamespace) -> Expression:
"""Return the pull `K / N psi(r) (x_j' - x_i')` of the velocity of the first member towards the second."""
first, second = namespace.first, namespace.second
separation = getattr(first, self.position) - getattr(second, self.position)
velocity = f'{self.position}_d1'
weight = namespace.K / self.n_agents * (1 + separation.dot_product(separation)) ** (-namespace.beta)
return weight * (getattr(second, velocity) - getattr(first, velocity))