Harmonic oscillator

A particle on a spring, optionally damped.

• Name: harmonic_oscillator
• Field: physics/mechanics
• Formalism: ordinary differential equations
• References: Hooke 1678

Factory

@register_model(
    field='physics/mechanics',
    tags=('mechanics', 'ode', 'oscillation'),
    references=('Hooke 1678',),
    defaults={'until': 20.0, 'dt': 0.05},
)
def harmonic_oscillator(
    x: float = 1.0,
    v: float = 0.0,
    *,
    k: float = 1.0,
    mass: float = 1.0,
    damping: float = 0.0,
) -> System:
    """Return a particle on a spring, optionally damped.

    Args:
        x:
            The initial displacement.

        v:
            The initial velocity.

        k:
            The stiffness.

        mass:
            The mass.

        damping:
            The linear damping coefficient.

    Returns:
        system:
            The model, with the particle `mass` in one dimension.

    Examples:
        >>> from skcomplex.models.physics.mechanics import harmonic_oscillator
        >>> harmonic_oscillator(k=4.0).interactions[0].k
        4.0
    """
    interactions: list[BaseInteraction] = [HarmonicPotentialInteraction('spring', target='mass', k=k)]
    if damping:
        interactions.append(DampingInteraction('damping', target='mass', gamma=damping))
    return System(
        elements=[MassPointElement('mass', x=x, x_d1=v, m=mass)],
        interactions=interactions,
        space=ContinuousSpace(n_dim=1),
    )

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

Equations

mass.x[0]: diff(x_mass_0(t), t, t) == (center_spring_0 - x_mass_0(t))*k_spring/m_mass