Nonlinear oscillators

Nonlinear oscillators: Duffing, Toda, Van der Pol, Stuart-Landau, Adler, impulses, billiards, and chaos.

AdlerPhaseElement(label, visual=None, *, phi=0.0, detuning=0.5, epsilon=1.0)

Bases: Element

Phase difference of a forced oscillator, phi' = detuning + epsilon sin(phi), the Adler equation.

The phase locks at a fixed point when |detuning| <= epsilon, and slips with the period 2 pi / sqrt(detuning^2 - epsilon^2) otherwise.

Parameters:

Name	Type	Description	Default
label	str	The label of the phase difference.	required
visual	VisualSpec | None	How the phase difference is drawn, or None for default_visual_.	None
phi	Annotated[float, StateSpec()]	The initial phase difference.	0.0
detuning	Annotated[float, ParameterSpec()]	The detuning of the natural frequency from the forcing.	0.5
epsilon	Annotated[float, ParameterSpec()]	The strength of the forcing.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the phase difference when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import AdlerPhaseElement
>>> AdlerPhaseElement('difference').detuning
0.5

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    phi: Annotated[float, StateSpec()] = 0.0,
    detuning: Annotated[float, ParameterSpec()] = 0.5,
    epsilon: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.phi = phi
    self.detuning = detuning
    self.epsilon = epsilon

DuffingElement(label, visual=None, *, x=1.0, x_d1=0.0, delta=0.3, alpha=-1.0, beta=1.0, gamma=0.5, omega=1.2)

Bases: Element

Displacement and parameters of the forced Duffing oscillator, whose acceleration a RateInteraction sets.

The acceleration is x'' = -delta x' - alpha x - beta x^3 + gamma cos(omega t).

For alpha < 0 < beta the potential has two wells, and a strong enough forcing drives the oscillator chaotically between them.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial displacement.	1.0
x_d1	float	The initial velocity.	0.0
delta	Annotated[float, ParameterSpec()]	The damping.	0.3
alpha	Annotated[float, ParameterSpec()]	The linear stiffness.	-1.0
beta	Annotated[float, ParameterSpec()]	The cubic stiffness.	1.0
gamma	Annotated[float, ParameterSpec()]	The amplitude of the forcing.	0.5
omega	Annotated[float, ParameterSpec()]	The frequency of the forcing.	1.2

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import DuffingElement
>>> DuffingElement('duffing').alpha
-1.0

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 1.0,
    x_d1: float = 0.0,
    delta: Annotated[float, ParameterSpec()] = 0.3,
    alpha: Annotated[float, ParameterSpec()] = -1.0,
    beta: Annotated[float, ParameterSpec()] = 1.0,
    gamma: Annotated[float, ParameterSpec()] = 0.5,
    omega: Annotated[float, ParameterSpec()] = 1.2,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.delta = delta
    self.alpha = alpha
    self.beta = beta
    self.gamma = gamma
    self.omega = omega

HenonHeilesElement(label, visual=VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None}), *, x=0.0, y=0.1, px=0.35, py=0.0)

Bases: Element

Star in the Hénon-Heiles potential V = (x^2 + y^2) / 2 + x^2 y - y^3 / 3, with the canonical momenta px, py.

The coordinates and the momenta are states of order 1, which a Hamiltonian interaction moves. Below the escape energy 1 / 6 the motion is regular for small energies and mostly chaotic near the escape energy.

Parameters:

Name	Type	Description	Default
label	str	The label of the star.	required
visual	VisualSpec | None	How the star is drawn.	VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None})
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.0
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.1
px	Annotated[float, StateSpec()]	The initial first momentum.	0.35
py	Annotated[float, StateSpec()]	The initial second momentum.	0.0

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import HenonHeilesElement
>>> HenonHeilesElement('star').px
0.35

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None}),
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.1,
    px: Annotated[float, StateSpec()] = 0.35,
    py: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.px = px
    self.py = py

ImpulsiveOscillatorElement(label, visual=None, *, x=0.0, x_d1=0.0, H=2.0, omega0=1.0)

Bases: Element

Displacement of an oscillator that rings freely between periodic impulses that kick its velocity.

A rate interaction gives the restoring force, and impulse interactions kick the velocity x_d1 by multiples of the amplitude H at their times, so the element holds the states and the parameters they read.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial displacement.	0.0
x_d1	float	The initial velocity.	0.0
H	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the impulses, half the jump of the velocity that one impulse of the symmetric train causes.	2.0
omega0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The natural frequency of a linear spring.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import ImpulsiveOscillatorElement
>>> ImpulsiveOscillatorElement('oscillator').H
2.0

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 0.0,
    x_d1: float = 0.0,
    H: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    omega0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.H = H
    self.omega0 = omega0

MushroomBilliardElement(label, visual=None, *, x=0.1, y=0.6, vx=0.6, vy=0.8, radius=1.0, width=0.4, height=1.0)

Bases: Element

Free particle in a mushroom table: a half disc cap of radius radius on a stem of half width width.

The particle flies straight, and reflects specularly off the cap, the floor of the cap beside the stem, the walls of the stem, and its floor at depth height. Its phase space mixes regular orbits that stay in the cap with chaotic orbits that enter the stem.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.1
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.6
vx	Annotated[float, StateSpec()]	The initial first velocity.	0.6
vy	Annotated[float, StateSpec()]	The initial second velocity.	0.8
radius	Annotated[float, ParameterSpec(domain=open(0, oo))]	The radius of the cap.	1.0
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The half width of the stem.	0.4
height	Annotated[float, ParameterSpec(domain=open(0, oo))]	The height of the stem.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import MushroomBilliardElement
>>> MushroomBilliardElement('ball').width
0.4

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_billiards.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.1,
    y: Annotated[float, StateSpec()] = 0.6,
    vx: Annotated[float, StateSpec()] = 0.6,
    vy: Annotated[float, StateSpec()] = 0.8,
    radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
    height: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.vx = vx
    self.vy = vy
    self.radius = radius
    self.width = width
    self.height = height

QuadrupoleElement(label, visual=VisualSpec('bodies', {'positions': ('q1', 'q2'), 'trail': None}), *, q1=0.0, q2=1.0, p1=0.5, p2=0.0, B=0.5)

Bases: Element

Surface of a nucleus in the quartic quadrupole potential, with the coordinates q1, q2 and momenta p1, p2.

The quadrupole deformation of the nuclear surface moves under the Hamiltonian H = (q² + p²) / 2 + B q1 (3 q2² - q1²) / √2 + q⁴ / 8, with q² = q1² + q2², which a Hamiltonian interaction turns into the equations of motion. The cubic term of strength B breaks the rotational symmetry and makes the motion chaotic at intermediate energies, while the quartic term keeps it bounded at every energy.

Parameters:

Name	Type	Description	Default
label	str	The label of the nucleus.	required
visual	VisualSpec | None	How the nucleus is drawn.	VisualSpec('bodies', {'positions': ('q1', 'q2'), 'trail': None})
q1	Annotated[float, StateSpec()]	The initial first coordinate.	0.0
q2	Annotated[float, StateSpec()]	The initial second coordinate.	1.0
p1	Annotated[float, StateSpec()]	The initial first momentum.	0.5
p2	Annotated[float, StateSpec()]	The initial second momentum.	0.0
B	Annotated[float, ParameterSpec()]	The strength of the cubic term.	0.5

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import QuadrupoleElement
>>> QuadrupoleElement('nucleus').B
0.5

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('bodies', {'positions': ('q1', 'q2'), 'trail': None}),
    *,
    q1: Annotated[float, StateSpec()] = 0.0,
    q2: Annotated[float, StateSpec()] = 1.0,
    p1: Annotated[float, StateSpec()] = 0.5,
    p2: Annotated[float, StateSpec()] = 0.0,
    B: Annotated[float, ParameterSpec()] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.q1 = q1
    self.q2 = q2
    self.p1 = p1
    self.p2 = p2
    self.B = B

SinaiBilliardElement(label, visual=VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'outline': _sinai_walls}), *, x=0.1, y=0.6, vx=0.6, vy=0.8, size=1.0, radius=0.4)

Bases: Element

Free particle in a square table of half side size with a central disc of radius radius, the Sinai billiard.

The particle flies straight, and its specular reflections off the four walls and the disc make its motion ergodic and chaotic.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the particle is drawn.	VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'outline': _sinai_walls})
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.1
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.6
vx	Annotated[float, StateSpec()]	The initial first velocity.	0.6
vy	Annotated[float, StateSpec()]	The initial second velocity.	0.8
size	Annotated[float, ParameterSpec(domain=open(0, oo))]	The half side of the square.	1.0
radius	Annotated[float, ParameterSpec(domain=open(0, oo))]	The radius of the disc.	0.4

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import SinaiBilliardElement
>>> SinaiBilliardElement('ball').radius
0.4

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_billiards.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec(
        'bodies',
        {'positions': ('x', 'y'), 'trail': None, 'outline': _sinai_walls},
    ),
    *,
    x: Annotated[float, StateSpec()] = 0.1,
    y: Annotated[float, StateSpec()] = 0.6,
    vx: Annotated[float, StateSpec()] = 0.6,
    vy: Annotated[float, StateSpec()] = 0.8,
    size: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.vx = vx
    self.vy = vy
    self.size = size
    self.radius = radius

ThreeDiskBilliardElement(label, visual=VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'outline': _three_disk_walls}), *, x=0.1, y=0.6, vx=0.6, vy=0.8, spacing=2.4, radius=0.4)

Bases: Element

Free particle scattered by three discs of radius radius at the corners of a triangle of side spacing.

The centers of the discs lie at the distance spacing / √3 from the origin, at the polar angles -π / 2, π / 6, and 5 π / 6. The particle flies straight and reflects specularly off the discs. The billiard is open, so almost every particle leaves the scattering region after a finite number of collisions, and the number of collisions depends on the incoming direction in a fractal way.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the particle is drawn.	VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'outline': _three_disk_walls})
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.1
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.6
vx	Annotated[float, StateSpec()]	The initial first velocity.	0.6
vy	Annotated[float, StateSpec()]	The initial second velocity.	0.8
spacing	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance between the centers of two discs.	2.4
radius	Annotated[float, ParameterSpec(domain=open(0, oo))]	The radius of the discs.	0.4

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import ThreeDiskBilliardElement
>>> ThreeDiskBilliardElement('ball').spacing
2.4

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_billiards.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec(
        'bodies',
        {'positions': ('x', 'y'), 'trail': None, 'outline': _three_disk_walls},
    ),
    *,
    x: Annotated[float, StateSpec()] = 0.1,
    y: Annotated[float, StateSpec()] = 0.6,
    vx: Annotated[float, StateSpec()] = 0.6,
    vy: Annotated[float, StateSpec()] = 0.8,
    spacing: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.4,
    radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.vx = vx
    self.vy = vy
    self.spacing = spacing
    self.radius = radius

TodaElement(label, visual=None, *, x=0.0, x_d1=0.0, d=0.1, a=0.88, omega=1.0)

Bases: Element

Displacement and parameters of the forced Toda oscillator, whose exponential spring is stiff on one side only.

The acceleration is x'' = -d x' - (e^x - 1) + a sin(omega t). The restoring force e^x - 1 grows exponentially for x > 0 and saturates at -1 for x < 0, so the oscillator is passive and asymmetric, and its resonance curve bends towards low frequencies.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial displacement.	0.0
x_d1	float	The initial velocity.	0.0
d	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The damping.	0.1
a	Annotated[float, ParameterSpec()]	The amplitude of the forcing.	0.88
omega	Annotated[float, ParameterSpec()]	The frequency of the forcing.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import TodaElement
>>> TodaElement('toda').d
0.1

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 0.0,
    x_d1: float = 0.0,
    d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    a: Annotated[float, ParameterSpec()] = 0.88,
    omega: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.d = d
    self.a = a
    self.omega = omega

VanDerPolElement(label, visual=None, *, x=2.0, x_d1=0.0, mu=1.0, A=0.0, omega=1.0)

Bases: Element

Van der Pol oscillator x'' = mu (1 - x^2) x' - x + A cos(omega t), with an optional periodic forcing.

The damping is negative for small amplitudes and positive for large ones, so the unforced oscillator, A = 0, settles on a limit cycle of amplitude near 2. A forcing of amplitude A and frequency omega locks it or drives it quasiperiodically or chaotically.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the oscillator is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial displacement.	2.0
x_d1	float	The initial velocity.	0.0
mu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the nonlinear damping.	1.0
A	Annotated[float, ParameterSpec()]	The amplitude of the forcing.	0.0
omega	Annotated[float, ParameterSpec()]	The frequency of the forcing.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the oscillator when it sets none.

Examples:

>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.physics.nonlinear_oscillators import VanDerPolElement
>>> system = System(elements=[VanDerPolElement('oscillator', x=0.1)])
>>> trajectory = Simulation(system, dt=0.1).prepare().run(until=50.0).trajectory_
>>> round(float(abs(trajectory.values('oscillator.x')[-300:]).max()), 1)
2.0

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 2.0,
    x_d1: float = 0.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    A: Annotated[float, ParameterSpec()] = 0.0,
    omega: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.mu = mu
    self.A = A
    self.omega = omega

adler_equation(phi=0.0, *, detuning=0.5, epsilon=1.0)

Return the Adler equation Δφ' = Δω + ε sin Δφ for the phase difference of an oscillator and its drive.

The phase difference locks at a fixed point when |Δω| < ε, inside the Arnold tongue. Outside it, the phase slips by 2π every 2π / sqrt(Δω² - ε²).

Parameters:

Name	Type	Description	Default
phi	float	The initial phase difference.	0.0
detuning	float	The frequency difference Δω.	0.5
epsilon	float	The strength of the coupling.	1.0

Returns:

Name	Type	Description
system	System	The system, with the phase difference difference.phi.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import adler_equation
>>> adler_equation(detuning=2.0).elements[0].detuning
2.0

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('synchronization', 'ode', 'phase', 'oscillation'),
    references=('Adler 1946', 'Datseris and Parlitz 2022'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def adler_equation(phi: float = 0.0, *, detuning: float = 0.5, epsilon: float = 1.0) -> System:
    """Return the Adler equation `Δφ' = Δω + ε sin Δφ` for the phase difference of an oscillator and its drive.

    The phase difference locks at a fixed point when `|Δω| < ε`, inside the Arnold tongue. Outside it, the phase slips
    by `2π` every `2π / sqrt(Δω² - ε²)`.

    Args:
        phi:
            The initial phase difference.

        detuning:
            The frequency difference `Δω`.

        epsilon:
            The strength of the coupling.

    Returns:
        system:
            The system, with the phase difference `difference.phi`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import adler_equation
        >>> adler_equation(detuning=2.0).elements[0].detuning
        2.0
    """
    return System(elements=[AdlerPhaseElement('difference', phi=phi, detuning=detuning, epsilon=epsilon)])

billiards(shape='sinai', x=0.1, y=0.6, vx=0.6, vy=0.8, *, size=1.0, disc_radius=0.4, cap_radius=1.0, width=0.4, height=1.0, spacing=2.4)

Return a particle moving freely on a billiard table and reflecting specularly off its walls.

The Sinai table is the square [-size, size]² with a disc of radius disc_radius at its center, and its motion is chaotic. The mushroom table is the upper half of a disc of radius cap_radius on a stem of half width width and height height, and its phase space mixes regular orbits that stay in the cap with chaotic orbits that enter the stem. The three-disk table is open: three discs of radius disc_radius sit at the corners of an equilateral triangle of side spacing centered at the origin, with one disc below the origin. A particle that comes from far away is scattered chaotically between the discs and leaves again, and the number of its collisions and its outgoing angle depend on its incoming direction through a fractal scattering function. Every bounce is an event, and the speed is conserved.

Parameters:

Name	Type	Description	Default
shape	str	'sinai', 'mushroom', or 'three_disk'.	'sinai'
x	float	The initial horizontal position.	0.1
y	float	The initial vertical position.	0.6
vx	float	The initial horizontal velocity.	0.6
vy	float	The initial vertical velocity.	0.8
size	float	The half side of the Sinai table.	1.0
disc_radius	float	The radius of the disc of the Sinai table and of the discs of the three-disk table.	0.4
cap_radius	float	The radius of the cap of the mushroom.	1.0
width	float	The half width of the stem of the mushroom.	0.4
height	float	The height of the stem of the mushroom.	1.0
spacing	float	The distance between the centers of the discs of the three-disk table.	2.4

Returns:

Name	Type	Description
system	System	The system, with the position particle.x, particle.y and the velocity particle.vx, particle.vy.

Raises:

Type	Description
ValidationError	The shape is unknown.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import billiards
>>> type(billiards('mushroom').elements[0]).__name__
'MushroomBilliardElement'

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_billiards.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('events', 'billiard', 'chaos', 'mixed phase space'),
    references=('Sinai 1970', 'Bunimovich 2001', 'Gaspard and Rice 1989', 'Datseris and Parlitz 2022'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def billiards(
    shape: str = 'sinai',
    x: float = 0.1,
    y: float = 0.6,
    vx: float = 0.6,
    vy: float = 0.8,
    *,
    size: float = 1.0,
    disc_radius: float = 0.4,
    cap_radius: float = 1.0,
    width: float = 0.4,
    height: float = 1.0,
    spacing: float = 2.4,
) -> System:
    """Return a particle moving freely on a billiard table and reflecting specularly off its walls.

    The Sinai table is the square `[-size, size]²` with a disc of radius `disc_radius` at its center, and its motion is
    chaotic. The mushroom table is the upper half of a disc of radius `cap_radius` on a stem of half width `width` and
    height `height`, and its phase space mixes regular orbits that stay in the cap with chaotic orbits that enter the
    stem. The three-disk table is open: three discs of radius `disc_radius` sit at the corners of an equilateral
    triangle of side `spacing` centered at the origin, with one disc below the origin. A particle that comes from far
    away is scattered chaotically between the discs and leaves again, and the number of its collisions and its outgoing
    angle depend on its incoming direction through a fractal scattering function. Every bounce is an event, and the
    speed is conserved.

    Args:
        shape:
            `'sinai'`, `'mushroom'`, or `'three_disk'`.

        x:
            The initial horizontal position.

        y:
            The initial vertical position.

        vx:
            The initial horizontal velocity.

        vy:
            The initial vertical velocity.

        size:
            The half side of the Sinai table.

        disc_radius:
            The radius of the disc of the Sinai table and of the discs of the three-disk table.

        cap_radius:
            The radius of the cap of the mushroom.

        width:
            The half width of the stem of the mushroom.

        height:
            The height of the stem of the mushroom.

        spacing:
            The distance between the centers of the discs of the three-disk table.

    Returns:
        system:
            The system, with the position `particle.x`, `particle.y` and the velocity `particle.vx`, `particle.vy`.

    Raises:
        ValidationError:
            The shape is unknown.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import billiards
        >>> type(billiards('mushroom').elements[0]).__name__
        'MushroomBilliardElement'
    """
    if shape == 'sinai':
        table = SinaiBilliardElement('particle', x=x, y=y, vx=vx, vy=vy, size=size, radius=disc_radius)
        walls = [
            _wall('right_wall', lambda s: s.target.x - s.target.size, 1),
            _wall('left_wall', lambda s: s.target.x + s.target.size, -1),
            _wall('top_wall', lambda s: s.target.y - s.target.size, 1),
            _wall('bottom_wall', lambda s: s.target.y + s.target.size, -1),
            _wall('disc', lambda s: s.target.x**2 + s.target.y**2 - s.target.radius**2, -1),
        ]
        return System(elements=[table], interactions=walls)
    if shape == 'mushroom':
        table = MushroomBilliardElement(
            'particle',
            x=x,
            y=y,
            vx=vx,
            vy=vy,
            radius=cap_radius,
            width=width,
            height=height,
        )
        walls = [
            _wall('cap', lambda s: s.target.x**2 + s.target.y**2 - s.target.radius**2, 1, lambda s: s.target.y > 0),
            _wall('cap_floor', lambda s: s.target.y, -1, lambda s: abs(s.target.x) > s.target.width),
            _wall('right_stem_wall', lambda s: s.target.x - s.target.width, 1, lambda s: s.target.y < 0),
            _wall('left_stem_wall', lambda s: s.target.x + s.target.width, -1, lambda s: s.target.y < 0),
            _wall('stem_floor', lambda s: s.target.y + s.target.height, -1),
        ]
        return System(elements=[table], interactions=walls)
    if shape == 'three_disk':
        table = ThreeDiskBilliardElement('particle', x=x, y=y, vx=vx, vy=vy, spacing=spacing, radius=disc_radius)
        walls = [_disc(f'disc_{index}', angle) for index, angle in enumerate(ANGLES)]
        return System(elements=[table], interactions=walls)
    error_msg = f'Unknown billiard shape {shape!r}. Use one of {", ".join(SHAPES)}.'
    raise ValidationError(error_msg)

duffing(x=1.0, v=0.0, *, delta=0.3, alpha=-1.0, beta=1.0, gamma=0.5, omega=1.2)

Return the forced Duffing oscillator x'' = -delta x' - alpha x - beta x^3 + gamma cos(omega t), chaotic here.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	1.0
v	float	The initial velocity.	0.0
delta	float	The damping.	0.3
alpha	float	The linear stiffness, negative for a double well.	-1.0
beta	float	The cubic stiffness.	1.0
gamma	float	The amplitude of the forcing.	0.5
omega	float	The angular frequency of the forcing.	1.2

Returns:

Name	Type	Description
system	System	The system, with the displacement duffing.x of order 2 and its velocity duffing.x_d1.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import duffing
>>> duffing(gamma=0.3).elements[0].gamma
0.3

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('chaos', 'ode', 'second order', 'forcing'),
    references=('Duffing 1918', 'Trzaska 2023'),
    defaults={'until': 200.0, 'dt': 0.05},
)
def duffing(
    x: float = 1.0,
    v: float = 0.0,
    *,
    delta: float = 0.3,
    alpha: float = -1.0,
    beta: float = 1.0,
    gamma: float = 0.5,
    omega: float = 1.2,
) -> System:
    """Return the forced Duffing oscillator `x'' = -delta x' - alpha x - beta x^3 + gamma cos(omega t)`, chaotic here.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        delta:
            The damping.

        alpha:
            The linear stiffness, negative for a double well.

        beta:
            The cubic stiffness.

        gamma:
            The amplitude of the forcing.

        omega:
            The angular frequency of the forcing.

    Returns:
        system:
            The system, with the displacement `duffing.x` of order 2 and its velocity `duffing.x_d1`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import duffing
        >>> duffing(gamma=0.3).elements[0].gamma
        0.3
    """
    acceleration = RateInteraction(
        'acceleration',
        target='duffing',
        state='x',
        rate=lambda s: (
            -s.target.delta * derivative(s.target.x, 1)
            - s.target.alpha * s.target.x
            - s.target.beta * s.target.x**3
            + s.target.gamma * cos(s.target.omega * s.t)
        ),
    )
    return System(
        elements=[
            DuffingElement('duffing', x=x, x_d1=v, delta=delta, alpha=alpha, beta=beta, gamma=gamma, omega=omega),
        ],
        interactions=[acceleration],
    )

forced_toda(x=0.0, v=0.0, *, d=0.1, a=0.88, omega=1.0)

Return the forced Toda oscillator x'' + d x' + e^x - 1 = a sin(omega t), a passive asymmetric oscillator.

The exponential spring is stiff for x > 0 and soft for x < 0, unlike the symmetric spring of the Duffing oscillator. Its free frequency falls from 1 as the amplitude grows, so the resonance curve bends to low driving frequencies and the response swings far to the soft side. Stronger driving gives coexisting responses and period doubling cascades to chaos as the driving period 2π / ω changes.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	0.0
v	float	The initial velocity.	0.0
d	float	The damping.	0.1
a	float	The amplitude of the forcing.	0.88
omega	float	The angular frequency of the forcing.	1.0

Returns:

Name	Type	Description
system	System	The system, with the displacement toda.x of order 2 and its velocity toda.x_d1.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import forced_toda
>>> forced_toda(a=1.2).elements[0].a
1.2

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('chaos', 'ode', 'second order', 'forcing', 'resonance'),
    references=('Toda 1967', 'Kurz and Lauterborn 1988', 'Datseris and Parlitz 2022'),
    defaults={'until': 300.0, 'dt': 0.05},
)
def forced_toda(
    x: float = 0.0,
    v: float = 0.0,
    *,
    d: float = 0.1,
    a: float = 0.88,
    omega: float = 1.0,
) -> System:
    """Return the forced Toda oscillator `x'' + d x' + e^x - 1 = a sin(omega t)`, a passive asymmetric oscillator.

    The exponential spring is stiff for `x > 0` and soft for `x < 0`, unlike the symmetric spring of the Duffing
    oscillator. Its free frequency falls from 1 as the amplitude grows, so the resonance curve bends to low driving
    frequencies and the response swings far to the soft side. Stronger driving gives coexisting responses and period
    doubling cascades to chaos as the driving period `2π / ω` changes.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        d:
            The damping.

        a:
            The amplitude of the forcing.

        omega:
            The angular frequency of the forcing.

    Returns:
        system:
            The system, with the displacement `toda.x` of order 2 and its velocity `toda.x_d1`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import forced_toda
        >>> forced_toda(a=1.2).elements[0].a
        1.2
    """
    return System(elements=[TodaElement('toda', x=x, x_d1=v, d=d, a=a, omega=omega)])

forced_van_der_pol(x=2.0, v=0.0, *, mu=1.0, A=1.0, omega=1.0)

Return the Van der Pol oscillator driven by a periodic force, x'' = μ (1 - x²) x' - x + A cos(ω t).

A strong enough force entrains the limit cycle, so the oscillator locks to the driving frequency inside Arnold tongues, and the motion can be quasiperiodic or chaotic outside them.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	2.0
v	float	The initial velocity.	0.0
mu	float	The strength of the nonlinear damping.	1.0
A	float	The amplitude of the force.	1.0
omega	float	The angular frequency of the force.	1.0

Returns:

Name	Type	Description
system	System	The system, with the displacement oscillator.x of order 2 and its velocity oscillator.x_d1.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import forced_van_der_pol
>>> forced_van_der_pol(A=0.5).elements[0].A
0.5

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('limit cycle', 'ode', 'second order', 'forcing', 'synchronization'),
    references=('van der Pol and van der Mark 1927', 'Trzaska 2023'),
    defaults={'until': 100.0, 'dt': 0.05},
)
def forced_van_der_pol(
    x: float = 2.0,
    v: float = 0.0,
    *,
    mu: float = 1.0,
    A: float = 1.0,
    omega: float = 1.0,
) -> System:
    """Return the Van der Pol oscillator driven by a periodic force, `x'' = μ (1 - x²) x' - x + A cos(ω t)`.

    A strong enough force entrains the limit cycle, so the oscillator locks to the driving frequency inside Arnold
    tongues, and the motion can be quasiperiodic or chaotic outside them.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        mu:
            The strength of the nonlinear damping.

        A:
            The amplitude of the force.

        omega:
            The angular frequency of the force.

    Returns:
        system:
            The system, with the displacement `oscillator.x` of order 2 and its velocity `oscillator.x_d1`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import forced_van_der_pol
        >>> forced_van_der_pol(A=0.5).elements[0].A
        0.5
    """
    return System(elements=[VanDerPolElement('oscillator', x=x, x_d1=v, mu=mu, A=A, omega=omega)])

henon_heiles(x=0.0, y=0.1, px=0.35, py=0.0)

Return the Hénon-Heiles system of a star in a galactic potential, whose motion is chaotic at higher energies.

The star moves under the Hamiltonian H = (px² + py²) / 2 + (x² + y²) / 2 + x² y - y³ / 3, which is conserved.

Parameters:

Name	Type	Description	Default
x	float	The initial first coordinate.	0.0
y	float	The initial second coordinate.	0.1
px	float	The initial first momentum.	0.35
py	float	The initial second momentum.	0.0

Returns:

Name	Type	Description
system	System	The system, with the coordinates star.x and star.y and their momenta star.px and star.py.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import henon_heiles
>>> henon_heiles(px=0.3).elements[0].px
0.3

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('chaos', 'ode', 'hamiltonian'),
    references=('Hénon and Heiles 1964', 'Datseris and Parlitz 2022'),
    defaults={'until': 200.0, 'dt': 0.05, 'solver': Yoshida4Solver(dt=0.01), 'invariants': {'star.energy': 1e-6}},
    observables=_get_henon_heiles_observables,
)
def henon_heiles(x: float = 0.0, y: float = 0.1, px: float = 0.35, py: float = 0.0) -> System:
    """Return the Hénon-Heiles system of a star in a galactic potential, whose motion is chaotic at higher energies.

    The star moves under the Hamiltonian `H = (px² + py²) / 2 + (x² + y²) / 2 + x² y - y³ / 3`, which is conserved.

    Args:
        x:
            The initial first coordinate.

        y:
            The initial second coordinate.

        px:
            The initial first momentum.

        py:
            The initial second momentum.

    Returns:
        system:
            The system, with the coordinates `star.x` and `star.y` and their momenta `star.px` and `star.py`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import henon_heiles
        >>> henon_heiles(px=0.3).elements[0].px
        0.3
    """
    return System(
        elements=[HenonHeilesElement('star', x=x, y=y, px=px, py=py)],
        interactions=[
            HamiltonianInteraction(
                'motion',
                target='star',
                coordinates=('x', 'y'),
                momenta=('px', 'py'),
                hamiltonian=_get_henon_heiles_hamiltonian,
            ),
        ],
    )

impulsive_oscillator(x=0.0, v=2.572, *, H=2.0, a=0.4165, gamma=0.0, spring='tangent', omega0=1.0)

Return an oscillator kicked by a periodic train of impulses of alternating sign.

The stiff spring gives x'' + tan(x) / cos²(x) = 2 H Σ_k (-1)^k δ(t - (2 k - 1) a), a strongly nonlinear oscillator whose free motion sin x = sin A sin(t / cos A) of amplitude A < π / 2 is known exactly. Between the impulses the oscillator rings freely, and every impulse is an event that kicks its velocity, by -2 H at the times (4 k + 1) a and by 2 H at the times (4 k + 3) a. A motion of amplitude A is periodic with the period 4 a of the impulses when an impulse reverses its velocity, which gives a = q arccos(H q² / sqrt((1 - q²) (1 - q² H²))) with q = cos A, for arctan H < A < π / 2. The defaults start the orbit of amplitude A = 1.2 for H = 2, from x = 0 with the speed tan A.

The temporal dipole shift gamma in (-1, 1) replaces the rectangular wave whose derivative gives the impulses with the derivative of an asymmetric saw-tooth wave: the negative impulses move to (4 k + 1 - gamma) a, the positive ones to (4 k + 3 + gamma) a, and both grow to 2 H / (1 - gamma²). The linear spring gives the harmonic oscillator x'' + omega0² x under the same impulses.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	0.0
v	float	The initial velocity.	2.572
H	float	The amplitude of the impulses.	2.0
a	float	A quarter of the period of the impulses.	0.4165
gamma	float	The temporal dipole shift of the impulses, in (-1, 1).	0.0
spring	str	The restoring force, 'tangent' for tan(x) / cos²(x) or 'linear' for omega0² x.	'tangent'
omega0	float	The natural frequency of the linear spring.	1.0

Returns:

Name	Type	Description
system	System	The system, with the displacement oscillator.x of order 2 and its velocity oscillator.x_d1, kicked by the impulses pull and push.

Raises:

Type	Description
ParameterError	The spring is unknown, the quarter period is not positive, or the shift is outside (-1, 1).

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import impulsive_oscillator
>>> impulsive_oscillator(a=0.5, gamma=0.5).interactions[1].start
0.25

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('ode', 'second order', 'impulses', 'events', 'forcing'),
    references=('Pilipchuk 2010', 'Trzaska 2023'),
    defaults={'until': 10.0, 'dt': 0.01, 'solver': ScipyODESolver(method='DOP853', rtol=1e-10, atol=1e-12)},
)
def impulsive_oscillator(
    x: float = 0.0,
    v: float = 2.572,
    *,
    H: float = 2.0,
    a: float = 0.4165,
    gamma: float = 0.0,
    spring: str = 'tangent',
    omega0: float = 1.0,
) -> System:
    """Return an oscillator kicked by a periodic train of impulses of alternating sign.

    The stiff spring gives `x'' + tan(x) / cos²(x) = 2 H Σ_k (-1)^k δ(t - (2 k - 1) a)`, a strongly nonlinear
    oscillator whose free motion `sin x = sin A sin(t / cos A)` of amplitude `A < π / 2` is known exactly. Between the
    impulses the oscillator rings freely, and every impulse is an event that kicks its velocity, by `-2 H` at the times
    `(4 k + 1) a` and by `2 H` at the times `(4 k + 3) a`. A motion of amplitude `A` is periodic with the period `4 a`
    of the impulses when an impulse reverses its velocity, which gives
    `a = q arccos(H q² / sqrt((1 - q²) (1 - q² H²)))` with `q = cos A`, for `arctan H < A < π / 2`. The defaults start
    the orbit of amplitude `A = 1.2` for `H = 2`, from `x = 0` with the speed `tan A`.

    The temporal dipole shift `gamma` in `(-1, 1)` replaces the rectangular wave whose derivative gives the impulses
    with the derivative of an asymmetric saw-tooth wave: the negative impulses move to `(4 k + 1 - gamma) a`, the
    positive ones to `(4 k + 3 + gamma) a`, and both grow to `2 H / (1 - gamma²)`. The linear spring gives the
    harmonic oscillator `x'' + omega0² x` under the same impulses.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        H:
            The amplitude of the impulses.

        a:
            A quarter of the period of the impulses.

        gamma:
            The temporal dipole shift of the impulses, in `(-1, 1)`.

        spring:
            The restoring force, `'tangent'` for `tan(x) / cos²(x)` or `'linear'` for `omega0² x`.

        omega0:
            The natural frequency of the linear spring.

    Returns:
        system:
            The system, with the displacement `oscillator.x` of order 2 and its velocity `oscillator.x_d1`, kicked by
            the impulses `pull` and `push`.

    Raises:
        ParameterError:
            The spring is unknown, the quarter period is not positive, or the shift is outside `(-1, 1)`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import impulsive_oscillator
        >>> impulsive_oscillator(a=0.5, gamma=0.5).interactions[1].start
        0.25
    """
    if spring not in SPRINGS:
        error_msg = f'Unknown spring {spring!r}. Use one of {", ".join(SPRINGS)}.'
        raise ParameterError(error_msg)
    if not a > 0 or not -1 < gamma < 1:
        error_msg = f'The impulses need a positive quarter period and a shift in (-1, 1), but got {a!r} and {gamma!r}.'
        raise ParameterError(error_msg)
    if spring == 'tangent':
        restoring = RateInteraction(
            'spring',
            target='oscillator',
            state='x',
            rate=lambda s: -tan(s.target.x) / cos(s.target.x) ** 2,
        )
    else:
        restoring = RateInteraction(
            'spring',
            target='oscillator',
            state='x',
            rate=lambda s: -s.target.omega0**2 * s.target.x,
        )
    jump = 2 / (1 - gamma**2)
    pull = {'x_d1': lambda s: s.target.x_d1 - jump * s.target.H}
    push = {'x_d1': lambda s: s.target.x_d1 + jump * s.target.H}
    return System(
        elements=[ImpulsiveOscillatorElement('oscillator', x=x, x_d1=v, H=H, omega0=omega0)],
        interactions=[
            restoring,
            ImpulseInteraction('pull', target='oscillator', reset=pull, period=4 * a, start=(1 - gamma) * a),
            ImpulseInteraction('push', target='oscillator', reset=push, period=4 * a, start=(3 + gamma) * a),
        ],
    )

nuclear_quadrupole(q1=0.0, q2=1.0, p1=0.5, p2=0.0, *, B=0.5)

Return the quadrupole vibrations of a nuclear surface, a Hamiltonian system of two degrees of freedom.

The deformation coordinates q = (q1, q2) and their momenta p move under the Hamiltonian H = (q² + p²) / 2 + B q1 (3 q2² - q1²) / √2 + q⁴ / 8, which is conserved. The cubic term has the threefold symmetry of the Hénon-Heiles potential, and the quartic term keeps every orbit bounded, so the motion is regular at low and high energies and largely chaotic in between.

Parameters:

Name	Type	Description	Default
q1	float	The initial first coordinate.	0.0
q2	float	The initial second coordinate.	1.0
p1	float	The initial first momentum.	0.5
p2	float	The initial second momentum.	0.0
B	float	The strength of the cubic term.	0.5

Returns:

Name	Type	Description
system	System	The system, with the coordinates nucleus.q1 and nucleus.q2, their momenta nucleus.p1 and nucleus.p2, and the observable nucleus.energy.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import nuclear_quadrupole
>>> nuclear_quadrupole(B=0.3).elements[0].B
0.3

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('chaos', 'ode', 'hamiltonian', 'nuclear physics'),
    references=('Baran and Raduta 1998', 'Datseris and Parlitz 2022'),
    defaults={'until': 200.0, 'dt': 0.05, 'solver': Yoshida4Solver(dt=0.01), 'invariants': {'nucleus.energy': 1e-6}},
    observables=_get_quadrupole_observables,
)
def nuclear_quadrupole(
    q1: float = 0.0,
    q2: float = 1.0,
    p1: float = 0.5,
    p2: float = 0.0,
    *,
    B: float = 0.5,
) -> System:
    """Return the quadrupole vibrations of a nuclear surface, a Hamiltonian system of two degrees of freedom.

    The deformation coordinates `q = (q1, q2)` and their momenta `p` move under the Hamiltonian
    `H = (q² + p²) / 2 + B q1 (3 q2² - q1²) / √2 + q⁴ / 8`, which is conserved. The cubic term has the threefold
    symmetry of the Hénon-Heiles potential, and the quartic term keeps every orbit bounded, so the motion is regular at
    low and high energies and largely chaotic in between.

    Args:
        q1:
            The initial first coordinate.

        q2:
            The initial second coordinate.

        p1:
            The initial first momentum.

        p2:
            The initial second momentum.

        B:
            The strength of the cubic term.

    Returns:
        system:
            The system, with the coordinates `nucleus.q1` and `nucleus.q2`, their momenta `nucleus.p1` and
            `nucleus.p2`, and the observable `nucleus.energy`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import nuclear_quadrupole
        >>> nuclear_quadrupole(B=0.3).elements[0].B
        0.3
    """
    return System(
        elements=[QuadrupoleElement('nucleus', q1=q1, q2=q2, p1=p1, p2=p2, B=B)],
        interactions=[
            HamiltonianInteraction(
                'motion',
                target='nucleus',
                coordinates=('q1', 'q2'),
                momenta=('p1', 'p2'),
                hamiltonian=_get_quadrupole_hamiltonian,
            ),
        ],
    )

stuart_landau(x=0.1, y=0.0, *, mu=1.0, omega=1.0, beta=0.0)

Return the Stuart-Landau oscillator z' = (μ + i ω) z - (1 + i β) |z|² z, the normal form of a Hopf bifurcation.

The complex amplitude z = x + i y is written as two real states. For μ > 0 every orbit except the origin settles on the circle of radius sqrt(μ), turning at the frequency ω - β μ.

Parameters:

Name	Type	Description	Default
x	float	The initial real part of the amplitude.	0.1
y	float	The initial imaginary part of the amplitude.	0.0
mu	float	The distance from the Hopf bifurcation.	1.0
omega	float	The linear frequency.	1.0
beta	float	The nonlinear frequency shift.	0.0

Returns:

Name	Type	Description
system	System	The system, with the real part oscillator.x and the imaginary part oscillator.y of the amplitude.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import stuart_landau
>>> stuart_landau(mu=0.5).elements[0].mu
0.5

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('limit cycle', 'ode', 'hopf bifurcation', 'oscillation'),
    references=('Gros 2015', 'Dilão 2023'),
    defaults={'until': 30.0, 'dt': 0.05},
)
def stuart_landau(
    x: float = 0.1,
    y: float = 0.0,
    *,
    mu: float = 1.0,
    omega: float = 1.0,
    beta: float = 0.0,
) -> System:
    """Return the Stuart-Landau oscillator `z' = (μ + i ω) z - (1 + i β) |z|² z`, the normal form of a Hopf bifurcation.

    The complex amplitude `z = x + i y` is written as two real states. For `μ > 0` every orbit except the origin
    settles on the circle of radius `sqrt(μ)`, turning at the frequency `ω - β μ`.

    Args:
        x:
            The initial real part of the amplitude.

        y:
            The initial imaginary part of the amplitude.

        mu:
            The distance from the Hopf bifurcation.

        omega:
            The linear frequency.

        beta:
            The nonlinear frequency shift.

    Returns:
        system:
            The system, with the real part `oscillator.x` and the imaginary part `oscillator.y` of the amplitude.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import stuart_landau
        >>> stuart_landau(mu=0.5).elements[0].mu
        0.5
    """
    return System(elements=[StuartLandauElement('oscillator', x=x, y=y, mu=mu, omega=omega, beta=beta)])

van_der_pol(x=2.0, v=0.0, *, mu=1.0)

Return the Van der Pol oscillator x'' = μ (1 - x²) x' - x, which settles on a limit cycle for μ > 0.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	2.0
v	float	The initial velocity.	0.0
mu	float	The strength of the nonlinear damping.	1.0

Returns:

Name	Type	Description
system	System	The system, with the displacement oscillator.x of order 2 and its velocity oscillator.x_d1.

Examples:

>>> from skcomplex.systems.physics.nonlinear_oscillators import van_der_pol
>>> van_der_pol(mu=3.0).elements[0].mu
3.0

Source code in src/skcomplex/systems/physics/nonlinear_oscillators/_oscillators.py

@register_system(
    field='physics/nonlinear_oscillators',
    tags=('limit cycle', 'ode', 'second order', 'oscillation'),
    references=('van der Pol 1926', 'Gros 2015'),
    defaults={'until': 50.0, 'dt': 0.05},
)
def van_der_pol(x: float = 2.0, v: float = 0.0, *, mu: float = 1.0) -> System:
    """Return the Van der Pol oscillator `x'' = μ (1 - x²) x' - x`, which settles on a limit cycle for `μ > 0`.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        mu:
            The strength of the nonlinear damping.

    Returns:
        system:
            The system, with the displacement `oscillator.x` of order 2 and its velocity `oscillator.x_d1`.

    Examples:
        >>> from skcomplex.systems.physics.nonlinear_oscillators import van_der_pol
        >>> van_der_pol(mu=3.0).elements[0].mu
        3.0
    """
    return System(elements=[VanDerPolElement('oscillator', x=x, x_d1=v, mu=mu)])