Driven damped oscillator
A mass on a spring with damping and a periodic force, m x'' = -k x - c x' + F cos(ω t).
- Name:
driven_damped_oscillator - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Hooke 1678
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'ode', 'second order', 'oscillation', 'forcing'),
references=('Hooke 1678',),
defaults={'until': 50.0, 'dt': 0.05},
)
def driven_damped_oscillator(
x: float = 1.0,
v: float = 0.0,
*,
m: float = 1.0,
k: float = 1.0,
c: float = 0.1,
F: float = 0.0,
omega: float = 1.0,
) -> System:
"""Return a mass on a spring with damping and a periodic force, `m x'' = -k x - c x' + F cos(ω t)`.
Args:
x:
The initial displacement.
v:
The initial velocity.
m:
The mass.
k:
The stiffness.
c:
The damping coefficient.
F:
The amplitude of the driving force.
omega:
The angular frequency of the driving force.
Returns:
system:
The model, with the mass `mass` in one dimension, whose displacement `mass.x` has order 2.
Examples:
>>> from skcomplex.models.physics.mechanics import driven_damped_oscillator
>>> driven_damped_oscillator(F=0.5).interactions[2].F
0.5
"""
return System(
elements=[MassPointElement('mass', x=x, x_d1=v, m=m)],
interactions=[
HarmonicPotentialInteraction('spring', target='mass', k=k),
DampingInteraction('damping', target='mass', gamma=c),
HarmonicForceInteraction('drive', target='mass', F=F, omega=omega),
],
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
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 HarmonicForceInteraction(RateInteraction):
"""Periodic force `F cos(omega t)` on a target, the second derivative `F cos(omega t) / m` of its position.
It is the rate `F cos(omega t)` with the mass as divisor, applied to every component of the position.
Args:
label:
The label of the interaction.
target:
The path of the element or group the force acts on.
F:
The amplitude of the force.
omega:
The angular frequency of the force.
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 skcomplex.models.physics.mechanics import HarmonicForceInteraction
>>> HarmonicForceInteraction('drive', target='mass', F=0.5, omega=2.0).omega
2.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
F: Annotated[float, ParameterSpec()] = 0.0,
omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
state: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, state=state, divisor=divisor)
self.F = F
self.omega = omega
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the force `F cos(omega t)`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The force, before the division by the mass.
"""
return namespace.F * cos(namespace.omega * namespace.t)
Equations
mass.x[0]: diff(x_mass_0(t), t, t) == (center_spring_0 - x_mass_0(t))*k_spring/m_mass + F_drive*cos(omega_drive*t)/m_mass - gamma_damping*diff(x_mass_0(t), t)/m_mass