Hasselmann stochastic climate
Hasselmann's stochastic climate, where fast weather noise drives a slow, damped temperature anomaly.
- Name:
hasselmann - Field:
physics/climate - Formalism: stochastic differential equations
- References: Hasselmann 1976; Ghil and Lucarini 2020
Factory
@register_model(
field='physics/climate',
tags=('sde', 'Ornstein-Uhlenbeck', 'red noise', 'stochastic climate'),
references=('Hasselmann 1976', 'Ghil and Lucarini 2020'),
defaults={'until': 100.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def hasselmann(T: float = 0.0, *, lam: float = 0.5, sigma: float = 0.3) -> System:
"""Return Hasselmann's stochastic climate, where fast weather noise drives a slow, damped temperature anomaly.
The anomaly follows `dT = -lam T dt + sigma dW`, an Ornstein-Uhlenbeck process whose red spectrum integrates white
weather noise, with the stationary variance `sigma^2 / (2 lam)`.
Args:
T:
The initial temperature anomaly.
lam:
The feedback rate that damps the anomaly.
sigma:
The amplitude of the weather noise.
Returns:
system:
The model, with the anomaly `ocean.T`.
Examples:
>>> from skcomplex.models.physics.climate import hasselmann
>>> hasselmann(lam=0.1).elements[0].lam
0.1
"""
feedback = RelaxationInteraction(
'feedback',
target='ocean',
state='T',
goal=0.0,
time=_get_damping_time,
noise=_get_weather,
)
return System(elements=[ClimateAnomalyElement('ocean', T=T, lam=lam, sigma=sigma)], interactions=[feedback])
Components
class ClimateAnomalyElement(Element):
"""Slow temperature anomaly of Hasselmann's stochastic climate, with its feedback rate and weather noise.
The relaxation of `hasselmann` damps the anomaly at the rate `lam` and drives it with weather noise of amplitude
`sigma`.
Args:
label:
The label of the anomaly.
visual:
How the element is drawn, or `None` for `default_visual_`.
T:
The initial temperature anomaly.
lam:
The feedback rate that damps the anomaly.
sigma:
The amplitude of the weather noise.
Examples:
>>> from skcomplex.models.physics.climate import ClimateAnomalyElement
>>> ClimateAnomalyElement('ocean').lam
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
T: Annotated[float, StateSpec()] = 0.0,
lam: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
) -> None:
super().__init__(label, visual=visual)
self.T = T
self.lam = lam
self.sigma = sigma
class RelaxationInteraction(RateInteraction):
"""Relaxation `(goal(ns) - x) / time(ns)` of a state `x` of a target towards a goal, with optional noise.
With a constant goal and time and an additive noise of amplitude `sigma`, the state follows the Ornstein-Uhlenbeck
process `dx = (goal - x) / time dt + sigma dW`, whose stationary variance is `sigma**2 time / 2`. Without noise it
is first-order smoothing, or adaptive expectation, towards the goal. A goal or a time function takes the namespace
of the interaction, as for `RateInteraction`, and a number stands for a constant. A goal that is a number applies
to every component of a vector state.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state relaxes.
state:
The name of the state of the target that relaxes.
goal:
A function of the namespace that returns the goal, or a number.
time:
A function of the namespace that returns the relaxation time, or a positive number.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
Raises:
ValidationError:
The state is not a name, or a function is neither callable nor a number.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RelaxationInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> smoothing = RelaxationInteraction('smoothing', target='belief', state='n', goal=10.0, time=0.5)
>>> system = System(elements=[CompartmentElement('belief', n=0.0)], interactions=[smoothing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('belief.n')[-1]), 3)
10.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
goal: Callable[[SimpleNamespace], object] | float = 0.0,
time: Callable[[SimpleNamespace], object] | float = 1.0,
noise: Callable[[SimpleNamespace], object] | float | None = None,
) -> None:
super().__init__(label, target=target, state=state, noise=noise)
self.goal = goal
self.time = time
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."""
super()._check_functions()
for name in ('goal', 'time'):
value = getattr(self, name)
if not (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)
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the goal, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.goal, namespace)
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.time, namespace)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `(goal - x) / time` of the state.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the state.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
goal = self._broadcast(self._get_goal(namespace), state)
return (goal - state) / self._get_time(namespace)
Equations
ocean.T: diff(T_ocean(t), t) == -lam_ocean*T_ocean(t)