Nonlinear oscillators

Nonlinear oscillators: Duffing, Van der Pol, Stuart-Landau, Adler, billiards, and Hamiltonian 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.models.physics.nonlinear_oscillators import AdlerPhaseElement
>>> AdlerPhaseElement('difference').detuning
0.5

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_adler_phase.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.models.physics.nonlinear_oscillators import DuffingElement
>>> DuffingElement('duffing').alpha
-1.0

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_duffing.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.models.physics.nonlinear_oscillators import HenonHeilesElement
>>> HenonHeilesElement('star').px
0.35

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_henon_heiles.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

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.models.physics.nonlinear_oscillators import MushroomBilliardElement
>>> MushroomBilliardElement('ball').width
0.4

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_mushroom_billiard.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

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.models.physics.nonlinear_oscillators import SinaiBilliardElement
>>> SinaiBilliardElement('ball').radius
0.4

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_sinai_billiard.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

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.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.nonlinear_oscillators import VanDerPolElement
>>> system = System(elements=[VanDerPolElement('oscillator', x=0.1)])
>>> trajectory = Simulation(system, dt=0.1).check().prepare().run(until=50.0).trajectory_
>>> round(float(abs(trajectory.values('oscillator.x')[-300:]).max()), 1)
2.0

Source code in src/skcomplex/models/physics/nonlinear_oscillators/_van_der_pol.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 model, with the phase difference difference.phi.

Examples:

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

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

@register_model(
    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 model, with the phase difference `difference.phi`.

    Examples:
        >>> from skcomplex.models.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)

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. Every bounce is an event, and the speed is conserved.

Parameters:

Name	Type	Description	Default
shape	str	'sinai' or 'mushroom'.	'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.	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

Returns:

Name	Type	Description
system	System	The model, 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.models.physics.nonlinear_oscillators import billiards
>>> type(billiards('mushroom').elements[0]).__name__
'MushroomBilliardElement'

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

@register_model(
    field='physics/nonlinear_oscillators',
    tags=('events', 'billiard', 'chaos', 'mixed phase space'),
    references=('Sinai 1970', 'Bunimovich 2001', '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,
) -> 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. Every bounce is an event, and the speed is conserved.

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

        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.

        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.

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

    Raises:
        ValidationError:
            The shape is unknown.

    Examples:
        >>> from skcomplex.models.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)
    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 model, with the displacement duffing.x of order 2 and its velocity duffing.x_d1.

Examples:

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

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

@register_model(
    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 model, with the displacement `duffing.x` of order 2 and its velocity `duffing.x_d1`.

    Examples:
        >>> from skcomplex.models.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_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 model, with the displacement oscillator.x of order 2 and its velocity oscillator.x_d1.

Examples:

>>> from skcomplex.models.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/models/physics/nonlinear_oscillators/_oscillators.py

@register_model(
    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 model, with the displacement `oscillator.x` of order 2 and its velocity `oscillator.x_d1`.

    Examples:
        >>> from skcomplex.models.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 model 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 model, with the coordinates star.x and star.y and their momenta star.px and star.py.

Examples:

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

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

@register_model(
    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 model 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 model, with the coordinates `star.x` and `star.y` and their momenta `star.px` and `star.py`.

    Examples:
        >>> from skcomplex.models.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,
            ),
        ],
    )

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 model, with the real part oscillator.x and the imaginary part oscillator.y of the amplitude.

Examples:

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

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

@register_model(
    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 model, with the real part `oscillator.x` and the imaginary part `oscillator.y` of the amplitude.

    Examples:
        >>> from skcomplex.models.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 model, with the displacement oscillator.x of order 2 and its velocity oscillator.x_d1.

Examples:

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

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

@register_model(
    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 model, with the displacement `oscillator.x` of order 2 and its velocity `oscillator.x_d1`.

    Examples:
        >>> from skcomplex.models.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)])