Nuclear quadrupole vibrations
The quadrupole vibrations of a nuclear surface, a Hamiltonian system of two degrees of freedom.
- Name:
nuclear_quadrupole - Field:
physics/nonlinear_oscillators - Formalism: ordinary differential equations
- References: Baran and Raduta 1998; Datseris and Parlitz 2022
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('nuclear_quadrupole')
entry = describe('nuclear_quadrupole')
trajectory = create_simulation('nuclear_quadrupole').prepare().run().trajectory_
Factory
@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,
),
],
)
Components
Elements:
QuadrupoleElement, fromElement
Interactions:
HamiltonianInteraction
class QuadrupoleElement(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.
Args:
label:
The label of the nucleus.
visual:
How the nucleus is drawn.
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.
Examples:
>>> from skcomplex.systems.physics.nonlinear_oscillators import QuadrupoleElement
>>> QuadrupoleElement('nucleus').B
0.5
"""
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
class HamiltonianInteraction(Interaction):
"""Hamilton's equations `q' = dH/dp` and `p' = -dH/dq` of the canonical pairs `(q, p)` of a Hamiltonian `H(ns)`.
The coordinates and the momenta are states of order 1, which may be angles, of every role that binds one element or
group, `target` by default, and the coordinates pair with the momenta in order. Every component of a vector state
is a coordinate or a momentum. The interaction returns one rate law per coordinate and per momentum, and the rates
of one state add, so the Hamiltonians of separate interactions sum: a coupling between two oscillators is an
interaction of its own whose Hamiltonian holds only the coupling energy. Dissipation or forcing is a rate
interaction on the momenta.
A Hamiltonian function takes the namespace of the interaction, with the namespaces of the roles, the time `t`, and
the symbols of the parameters a subclass declares.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states are the canonical pairs.
coordinates:
The name of the state of order 1 that is a coordinate, or the names of several states.
momenta:
The name of the state of order 1 that is the momentum of the coordinate, or the names of several states in
the order of the coordinates.
hamiltonian:
A function of the namespace that returns the Hamiltonian.
Raises:
ValidationError:
A state is not a name, the coordinates and the momenta differ in count, the Hamiltonian is not callable,
or a coordinate and its momentum differ in size or are not of order 1.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.interactions import HamiltonianInteraction
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import StateSpec
>>> class OscillatorElement(Element):
... def __init__(self, label, q: Annotated[float, StateSpec()] = 1.0, p: Annotated[float, StateSpec()] = 0):
... super().__init__(label)
... self.q = q
... self.p = p
>>> energy = lambda s: (s.target.p**2 + s.target.q**2) / 2
>>> spring = HamiltonianInteraction('spring', target='mass', hamiltonian=energy)
>>> system = System(elements=[OscillatorElement('mass')], interactions=[spring])
>>> trajectory = Simulation(system, dt=0.01).prepare().run(until=3.14159).trajectory_
>>> round(float(trajectory.values('mass.q')[-1]), 3)
-1.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('coordinates', 'momenta'))],
coordinates: str | Sequence[str] = 'q',
momenta: str | Sequence[str] = 'p',
hamiltonian: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.target = target
self.coordinates = coordinates
self.momenta = momenta
self.hamiltonian = hamiltonian
def _get_names(self, names: str | Sequence[str]) -> tuple[str, ...]:
"""Return the names of the states of one side of the canonical pairs."""
return (names,) if isinstance(names, str) else tuple(names)
def _check(self) -> None:
"""Raise an error for a state that is not a name, an unpaired state, or a Hamiltonian that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, the coordinates and the momenta differ in count, or
the Hamiltonian is not callable.
"""
super()._check()
for side in ('coordinates', 'momenta'):
names = self._get_names(getattr(self, side))
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'The {side} of {self.label!r} are {getattr(self, side)!r}. Give names of states.'
raise ValidationError(error_msg)
coordinates, momenta = self._get_names(self.coordinates), self._get_names(self.momenta)
if len(coordinates) != len(momenta):
unpaired = coordinates[len(momenta) :] or momenta[len(coordinates) :]
error_msg = (
f'{type(self).__name__} {self.label!r} pairs the coordinates {coordinates} with the momenta {momenta}, '
f'so {", ".join(unpaired)} has no partner. Give one momentum per coordinate.'
)
raise ValidationError(error_msg)
if self.hamiltonian is not None and not callable(self.hamiltonian):
error_msg = f'The hamiltonian of {self.label!r} is {self.hamiltonian!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_hamiltonian(self, namespace: SimpleNamespace) -> object:
"""Return the Hamiltonian, the hook a subclass overrides to derive it from its parameters."""
if self.hamiltonian is None:
error_msg = f'{type(self).__name__} {self.label!r} has no Hamiltonian. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.hamiltonian(namespace)
def _get_members(self, namespace: SimpleNamespace) -> list[SimpleNamespace]:
"""Return the namespaces of the members whose states are canonical pairs, those of every role of one path."""
return [
getattr(namespace, name)
for name, spec in self._get_specs().items()
if isinstance(spec, RoleSpec) and not spec.variadic and not spec.noise
]
def _get_state(self, member: SimpleNamespace, name: str) -> Target:
"""Return a canonical state of a member, after checking that it has order 1."""
if hasattr(member, f'{name}_d1'):
error_msg = (
f'{type(self).__name__} {self.label!r} reads the state {name} of order above 1, but canonical '
'coordinates and momenta have order 1. Declare it with `StateSpec(order=1)`, or derive the equations '
'of a coordinate of order 2 with `EulerLagrangeInteraction`.'
)
raise ValidationError(error_msg)
state: Target = getattr(member, name)
return state
def _get_pairs(self, namespace: SimpleNamespace) -> list[tuple[Target, Target]]:
"""Return the canonical pairs of every member, after checking that each coordinate matches its momentum."""
pairs = []
names = list(zip(self._get_names(self.coordinates), self._get_names(self.momenta), strict=True))
for member in self._get_members(namespace):
for coordinate, momentum in names:
q, p = self._get_state(member, coordinate), self._get_state(member, momentum)
if len(self._get_components(q)) != len(self._get_components(p)):
error_msg = (
f'The coordinate {coordinate} and the momentum {momentum} of {self.label!r} differ in size. '
'Give a momentum of the size of its coordinate.'
)
raise ValidationError(error_msg)
pairs.append((q, p))
return pairs
def _get_components(self, state: Target) -> list[Expression]:
"""Return the components of a state, the state itself for a scalar."""
return list(state) if isinstance(state, FreeModuleElement) else [state]
def _get_derivative(self, hamiltonian: Expression, state: Target) -> Target:
"""Return the derivative of the Hamiltonian with respect to a state, every component of a vector state."""
if isinstance(state, FreeModuleElement):
return vector(SR, [hamiltonian.diff(component) for component in state])
return hamiltonian.diff(state)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate of every coordinate, `dH/dp`, and of every momentum, `-dH/dq`.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The rate laws of the coordinates and the momenta.
"""
pairs = self._get_pairs(namespace)
hamiltonian = SR(self._get_hamiltonian(namespace))
laws: list[BasePart] = []
for q, p in pairs:
laws.append(RateLaw(q, self._get_derivative(hamiltonian, p)))
laws.append(RateLaw(p, -self._get_derivative(hamiltonian, q)))
return laws
Equations
nucleus.q1: diff(q1_nucleus(t), t) == p1_nucleus(t)
nucleus.q2: diff(q2_nucleus(t), t) == p2_nucleus(t)
nucleus.p1: diff(p1_nucleus(t), t) == sqrt(2)*B_nucleus*q1_nucleus(t)^2 + 1/2*sqrt(2)*(q1_nucleus(t)^2 - 3*q2_nucleus(t)^2)*B_nucleus - 1/2*(q1_nucleus(t)^2 + q2_nucleus(t)^2)*q1_nucleus(t) - q1_nucleus(t)
nucleus.p2: diff(p2_nucleus(t), t) == -3*sqrt(2)*B_nucleus*q1_nucleus(t)*q2_nucleus(t) - 1/2*(q1_nucleus(t)^2 + q2_nucleus(t)^2)*q2_nucleus(t) - q2_nucleus(t)