Strongly nonlinear oscillator under periodic impulses
An oscillator kicked by a periodic train of impulses of alternating sign.
- Name:
impulsive_oscillator - Field:
physics/nonlinear_oscillators - Formalism: ordinary differential equations
- References: Pilipchuk 2010; Trzaska 2023
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('impulsive_oscillator')
entry = describe('impulsive_oscillator')
trajectory = create_simulation('impulsive_oscillator').prepare().run().trajectory_
Factory
@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),
],
)
Components
Elements:
ImpulsiveOscillatorElement, fromElement
Interactions:
RateInteractionImpulseInteraction
class ImpulsiveOscillatorElement(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.
Args:
label:
The label of the oscillator.
visual:
How the oscillator is drawn, or `None` for `default_visual_`.
x:
The initial displacement.
x_d1:
The initial velocity.
H:
The amplitude of the impulses, half the jump of the velocity that one impulse of the symmetric train
causes.
omega0:
The natural frequency of a linear spring.
Attributes:
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
"""
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
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.components.elements import CompartmentElement
>>> from skcomplex.core.components.interactions import RateInteraction
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw, NoiseLaw, FractionalNoiseLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'divisor'))],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
class ImpulseInteraction(Interaction):
"""Impulses that set states of the target to new values at given times, or at periodic times.
The impulses fall at every time of `times` and, with a `period`, at `start + j * period` for `j` from 0, `count`
times or without end. An impulse at the initial time of a run does not fire. The new values are functions of the
namespace, which holds `namespace.target`, the states and parameters of the target, and the time `namespace.t`. An
impulse sets the states of every member of a group at once, and the trajectory records every impulse among its
events, under the label of the interaction.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the impulses set.
reset:
A mapping from the name of a state of the target to the function `f(namespace)` of its new value.
times:
The times of the impulses.
period:
The period of the periodic impulses, or `None` for none.
start:
The time of the first periodic impulse.
count:
The number of periodic impulses, or `None` for impulses without end.
Raises:
ValidationError:
The impulses set no state, a new value is not a function, the impulses have no time, the period is not
positive, or the count is not a positive integer.
Examples:
>>> from skcomplex.core.components.interactions import ImpulseInteraction
>>> reset = {'c': lambda s: s.target.c + 1.0}
>>> dose = ImpulseInteraction('dose', target='body', reset=reset, period=8.0, count=3)
>>> dose.start
0.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (CrossingLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('reset',))],
reset: Mapping[str, Callable[[SimpleNamespace], object]],
times: Sequence[float] = (),
period: float | None = None,
start: float = 0.0,
count: int | None = None,
) -> None:
super().__init__(label)
self.target = target
self.reset = reset
self.times = times
self.period = period
self.start = start
self.count = count
def _check(self) -> None:
"""Raise an error for impulses without states or times, or with a period or a count they cannot take.
Raises:
ValidationError:
The impulses set no state, a new value is not a function, the impulses have no time, the period is not
positive, or the count is not a positive integer.
"""
super()._check()
if not self.reset or not all(callable(function) for function in self.reset.values()):
error_msg = (
f'The impulses {self.label!r} set {self.reset!r}. Map state names to functions of the namespace.'
)
raise ValidationError(error_msg)
if not self.times and self.period is None:
error_msg = f'The impulses {self.label!r} have no time. Give `times`, or a `period`.'
raise ValidationError(error_msg)
if self.period is not None and not self.period > 0:
error_msg = f'The period of the impulses {self.label!r} is {self.period!r}. Give a positive number.'
raise ValidationError(error_msg)
if self.count is not None and (
isinstance(self.count, bool) or not isinstance(self.count, int) or self.count < 1
):
error_msg = f'The count of the impulses {self.label!r} is {self.count!r}. Give a positive integer, or None.'
raise ValidationError(error_msg)
def _get_periodic_guard(self, time: Expression, period: float) -> Expression:
"""Return a guard that crosses zero at every periodic impulse and stays away from zero outside them."""
guard = where(time < self.start - period / 2, -1, sin(pi * (time - self.start) / period))
if self.count is None:
return guard
last = self.count - SR(1) / 2
return where(time > self.start + last * period, sin(pi * last), guard)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing laws of time at the impulses.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target` and the time `namespace.t`.
Returns:
laws:
One crossing law per given time and one for the periodic impulses.
"""
reset = {
getattr(self._get_member(namespace, 'target'), name): SR(f(namespace)) for name, f in self.reset.items()
}
laws: list[BaseLaw] = [CrossingLaw(namespace.t - time, reset, 1, label=self.label) for time in self.times]
if self.period is not None:
laws.append(CrossingLaw(self._get_periodic_guard(namespace.t, self.period), reset, 0, label=self.label))
return laws
Equations
oscillator.x: diff(x_oscillator(t), t, t) == -tan(x_oscillator(t))/cos(x_oscillator(t))^2