Magnetic pendulum
A magnetic bob on a string swinging above magnets on a circle, which comes to rest at one of them.
- Name:
magnetic_pendulum - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Datseris and Parlitz 2022
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'ode', 'multistability', 'fractal basins'),
references=('Datseris and Parlitz 2022',),
defaults={'until': 50.0, 'dt': 0.05},
)
def magnetic_pendulum(
x: float = 1.0,
y: float = 1.5,
*,
n_magnets: int = 3,
omega: float = 1.0,
q: float = 0.2,
gamma: float = 1.0,
d: float = 0.3,
) -> System:
"""Return a magnetic bob on a string swinging above magnets on a circle, which comes to rest at one of them.
The bob follows `x'' = -ω² x - q x' - Σ gamma (x - x_i) / d_i³` with `d_i = sqrt(|x - x_i|² + d²)`. The basins of
the magnets have fractal boundaries.
Args:
x:
The initial horizontal position of the bob.
y:
The initial vertical position of the bob.
n_magnets:
The number of magnets, evenly spaced on the unit circle.
omega:
The natural frequency of the pendulum.
q:
The friction coefficient.
gamma:
The strength of every magnet.
d:
The height of the bob above the plane of the magnets.
Returns:
system:
The model, with the bob `bob` and its position `bob.x` in two dimensions.
Examples:
>>> from skcomplex.models.physics.mechanics import magnetic_pendulum
>>> len(magnetic_pendulum(n_magnets=4).interactions)
6
"""
angles = 2 * np.pi * np.arange(n_magnets) / n_magnets
magnets = [
MagnetInteraction(
f'magnet_{index}',
target='bob',
center=[round(float(np.cos(angle)), 12), round(float(np.sin(angle)), 12)],
gamma=gamma,
d=d,
)
for index, angle in enumerate(angles)
]
return System(
elements=[MassPointElement('bob', x=[x, y], x_d1=[0.0, 0.0])],
interactions=[
HarmonicPotentialInteraction('string', target='bob', k=omega**2),
DampingInteraction('friction', target='bob', gamma=q),
*magnets,
],
space=ContinuousSpace(n_dim=2),
)
Components
class MassPointElement(PointElement):
"""Point mass, a point whose position `x` has order 2, with its velocity `x_d1` and a mass `m`.
The forces on the point divide by the mass through their divisor `'m'`.
Args:
label:
The label of the point.
x:
The initial position. A number applies to every component.
x_d1:
The initial velocity. A number applies to every component.
visual:
How the point is drawn, or `None` for `default_visual_`.
m:
The mass.
Attributes:
default_visual_:
Points at the positions.
Examples:
>>> from skcomplex.models.physics.common.particles import MassPointElement
>>> MassPointElement('ball', x=[0.0, 1.0], m=2.0).m
2.0
"""
def __init__(
self,
label: str,
x: npt.ArrayLike = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
*,
m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, x=x, x_d1=x_d1, visual=visual)
self.m = m
class HarmonicPotentialInteraction(PotentialEnergyInteraction):
"""Harmonic potential that pulls a particle towards a center.
The energy is `k * |x - center|**2 / 2`.
Args:
label:
The label of the interaction.
target:
The path of the particle.
k:
The stiffness.
center:
The center. A number applies to every component.
states:
The name of the vector state that holds the position of the particle.
divisor:
The name of the state or parameter that holds the mass of the particle.
Examples:
>>> from skcomplex.models.physics.common.particles import HarmonicPotentialInteraction
>>> HarmonicPotentialInteraction('trap', target='ball', k=2.0).k
2.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
center: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
states: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, states=states, divisor=divisor)
self.k = k
self.center = center
def _get_center(self, s: SimpleNamespace) -> object:
"""Return the center, the hook a subclass overrides to move it."""
return s.center
def _get_potential(self, s: SimpleNamespace) -> Expression:
"""Return `k * |x - center|**2 / 2`.
Args:
s:
The namespace of the symbols the energy can use.
Returns:
energy:
The potential energy as a symbolic expression.
"""
displacement = getattr(s.target, self._get_names()[0]) - self._get_center(s)
return s.k * displacement.dot_product(displacement) / 2
class DampingInteraction(RateInteraction):
"""Linear damping, a force `-gamma * v` proportional and opposite to the velocity `v` of a target.
The velocity is the first derivative of the position, a state of order 2, whose second derivative changes by
`-gamma * v / m` with the mass as divisor. The LAMMPS backend runs it as `fix viscous`.
Args:
label:
The label of the interaction.
target:
The path of the element or group whose velocity is damped.
gamma:
The damping coefficient.
state:
The name of the state of order 2 that holds the position of the target.
divisor:
The name of the state or parameter that holds the mass of the target.
Examples:
>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.models.physics.common.particles import DampingInteraction
>>> target = SimpleNamespace(x=var('x'), m=2)
>>> damping = DampingInteraction('air', target='ball', gamma=0.1).check().prepare()
>>> [law.value for law in damping.run(SimpleNamespace(target=target, gamma=var('gamma'))).laws_]
[-1/2*gamma*derivative(x, 1)]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
state: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, state=state, divisor=divisor)
self.gamma = gamma
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the drag `-gamma * v`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The drag, before the division by the mass.
"""
return -namespace.gamma * velocity_of(getattr(namespace.target, self.state))
class MagnetInteraction(PotentialEnergyInteraction):
"""Magnet below the plane of a pendulum bob, which attracts it with energy `-gamma / sqrt(|x - center|**2 + d**2)`.
Args:
label:
The label of the interaction.
target:
The path of the bob.
center:
The position of the magnet in the plane.
gamma:
The strength of the magnet.
d:
The height of the bob above the plane of the magnets.
Examples:
>>> from skcomplex.models.physics.mechanics import MagnetInteraction
>>> MagnetInteraction('magnet', target='bob', center=[1.0, 0.0]).gamma
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
center: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
d: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
) -> None:
super().__init__(label, target=target)
self.center = center
self.gamma = gamma
self.d = d
def _get_potential(self, s: SimpleNamespace) -> Expression:
"""Return `-gamma / sqrt(|x - center|**2 + d**2)`.
Args:
s:
The namespace of the symbols the energy can use.
Returns:
energy:
The potential energy as a symbolic expression.
"""
displacement = s.target.x - s.center
return -s.gamma / sqrt(displacement.dot_product(displacement) + s.d**2)
Equations
bob.x[0]: diff(x_bob_0(t), t, t) == (center_string_0 - x_bob_0(t))*k_string/m_bob - gamma_friction*diff(x_bob_0(t), t)/m_bob + (center_magnet_0_0 - x_bob_0(t))*gamma_magnet_0/(((center_magnet_0_0 - x_bob_0(t))^2 + (center_magnet_0_1 - x_bob_1(t))^2 + d_magnet_0^2)^(3/2)*m_bob) + (center_magnet_1_0 - x_bob_0(t))*gamma_magnet_1/(((center_magnet_1_0 - x_bob_0(t))^2 + (center_magnet_1_1 - x_bob_1(t))^2 + d_magnet_1^2)^(3/2)*m_bob) + (center_magnet_2_0 - x_bob_0(t))*gamma_magnet_2/(((center_magnet_2_0 - x_bob_0(t))^2 + (center_magnet_2_1 - x_bob_1(t))^2 + d_magnet_2^2)^(3/2)*m_bob)
bob.x[1]: diff(x_bob_1(t), t, t) == (center_string_1 - x_bob_1(t))*k_string/m_bob - gamma_friction*diff(x_bob_1(t), t)/m_bob + (center_magnet_0_1 - x_bob_1(t))*gamma_magnet_0/(((center_magnet_0_0 - x_bob_0(t))^2 + (center_magnet_0_1 - x_bob_1(t))^2 + d_magnet_0^2)^(3/2)*m_bob) + (center_magnet_1_1 - x_bob_1(t))*gamma_magnet_1/(((center_magnet_1_0 - x_bob_0(t))^2 + (center_magnet_1_1 - x_bob_1(t))^2 + d_magnet_1^2)^(3/2)*m_bob) + (center_magnet_2_1 - x_bob_1(t))*gamma_magnet_2/(((center_magnet_2_0 - x_bob_0(t))^2 + (center_magnet_2_1 - x_bob_1(t))^2 + d_magnet_2^2)^(3/2)*m_bob)