|
| 1 | +""" |
| 2 | +This file contains numerical solvers for the Kuramoto model. You can add your own solvers here (write your own decorator |
| 3 | +and add it to include.solvers if you wish). |
| 4 | +""" |
| 5 | +from typing import Union |
| 6 | +import torch |
| 7 | + |
| 8 | +# Import the solver module (located in `include`) |
| 9 | +import sys |
| 10 | +from os.path import dirname as up |
| 11 | +from dantro._import_tools import import_module_from_path |
| 12 | +sys.path.append(up(up(__file__))) |
| 13 | +include = import_module_from_path(mod_path=up(up(__file__)), mod_str="include") |
| 14 | +from include.solvers import torchdiffeq_solver |
| 15 | + |
| 16 | + |
| 17 | +def Kuramoto_rhs( |
| 18 | + t, |
| 19 | + state, |
| 20 | + *, |
| 21 | + adjacency_matrix: torch.Tensor, |
| 22 | + eigen_frequencies: torch.Tensor, |
| 23 | + kappa: Union[torch.Tensor, float], |
| 24 | + beta: Union[torch.Tensor, float], |
| 25 | + alpha: Union[torch.Tensor, float] = 0.0, |
| 26 | +): |
| 27 | + """ |
| 28 | + Compute the right-hand side (RHS) of the Kuramoto model. |
| 29 | +
|
| 30 | + Parameters |
| 31 | + ---------- |
| 32 | + t : float or torch.Tensor |
| 33 | + Current time (unused in the autonomous Kuramoto equations, but kept |
| 34 | + for compatibility with generic ODE solver decorators). |
| 35 | + state : torch.Tensor |
| 36 | + Current system state. |
| 37 | + - If `alpha == 0` (first-order dynamics): shape (N, 1), |
| 38 | + containing oscillator phases θ_i. |
| 39 | + - If `alpha != 0` (second-order dynamics): shape (N, 2), |
| 40 | + with [:, 0] = phases θ_i and [:, 1] = velocities dθ_i/dt. |
| 41 | + adjacency_matrix : torch.Tensor, shape (N, N) |
| 42 | + Coupling matrix describing network connections. |
| 43 | + eigen_frequencies : torch.Tensor, shape (N, 1) |
| 44 | + Natural frequencies ω_i of the oscillators. |
| 45 | + kappa : float |
| 46 | + Global coupling strength. |
| 47 | + alpha : float |
| 48 | + Inertia parameter. If 0, reduces to the standard first-order |
| 49 | + Kuramoto model. Otherwise, adds a second-order (inertial) term. |
| 50 | + beta : float |
| 51 | + Damping parameter. Used in both first- and second-order cases. |
| 52 | + sigma : float |
| 53 | + Noise strength (unused here — stochastic increments should be |
| 54 | + added at the solver level, not inside the deterministic RHS). |
| 55 | + device : torch.device |
| 56 | + Torch device for tensor operations. |
| 57 | +
|
| 58 | + Returns |
| 59 | + ------- |
| 60 | + torch.Tensor |
| 61 | + Derivatives of the state with the same shape as `state`: |
| 62 | + - First-order (alpha == 0): shape (N, 1), dθ/dt |
| 63 | + - Second-order (alpha != 0): shape (N, 2), |
| 64 | + [dθ/dt, d²θ/dt²] |
| 65 | + """ |
| 66 | + phases = state[:, 0] |
| 67 | + if alpha != 0: # second-order case |
| 68 | + velocities = state[:, 1] |
| 69 | + else: |
| 70 | + velocities = None |
| 71 | + |
| 72 | + # Pairwise phase differences |
| 73 | + diffs = torch.sin(phases - phases.reshape((len(phases),))) |
| 74 | + |
| 75 | + # Coupling contribution |
| 76 | + coupling = torch.matmul(kappa * adjacency_matrix, diffs).diag() |
| 77 | + |
| 78 | + if alpha == 0: |
| 79 | + # First-order Kuramoto |
| 80 | + dtheta = (eigen_frequencies.squeeze() + coupling) / beta |
| 81 | + return torch.stack([dtheta], dim=1) |
| 82 | + else: |
| 83 | + # Second-order Kuramoto |
| 84 | + dtheta = velocities |
| 85 | + dvel = (eigen_frequencies.squeeze() + coupling - beta * velocities) / alpha |
| 86 | + return torch.stack([dtheta, dvel], dim=1) |
| 87 | + |
| 88 | +# Euler solver |
| 89 | +@torchdiffeq_solver(method="euler", adjoint=False) |
| 90 | +def Kuramoto_euler(t, state, **params): |
| 91 | + return Kuramoto_rhs(t, state, **params) |
| 92 | + |
| 93 | + |
| 94 | +# Dopri5 solver |
| 95 | +@torchdiffeq_solver(method="dopri5", adjoint=False) |
| 96 | +def Kuramoto_dopri5(t, state, **params): |
| 97 | + return Kuramoto_rhs(t, state, **params) |
| 98 | + |
| 99 | + |
| 100 | +# Runge-Kutta 4th order solver |
| 101 | +@torchdiffeq_solver(method="rk4", adjoint=False) |
| 102 | +def Kuramoto_rk4_adj(t, state, **params): |
| 103 | + return Kuramoto_rhs(t, state, **params) |
| 104 | + |
0 commit comments