catlog/simulate/ode/
kuramoto.rs

1//! Kuramoto model of synchronization of coupled oscillators.
2//!
3//! The first-order and second-order Kuramoto models are both described in
4//! Section 2.1 of [Nitzbon et al 2017](crate::refs::NitzbonNetworkStability),
5//! among other sources.
6
7use nalgebra::{DMatrix, DVector};
8
9use super::ODESystem;
10
11/// Differential order of a Kuramoto system.
12#[derive(Clone, Copy, Debug, PartialEq, Eq)]
13pub enum KuramotoOrder {
14    /// First-order Kuramoto system.
15    First,
16
17    /// Second-order Kuramoto system.
18    ///
19    /// Reduced to a first-order ODE system by introducing extra state variables
20    /// for the angular frequencies `ω_i`.
21    Second,
22}
23
24/// Kuramoto system of ODEs.
25///
26/// The state variables of the system are the phase angles `ϕ_i` and, in the
27/// second-order case, also the angular frequencies `ω_i`. The notation here
28/// follows [Nitzbon et al 2017](crate::refs::NitzbonNetworkStability).
29#[derive(Clone, Debug, PartialEq)]
30pub struct KuramotoSystem {
31    /// Differential order of the system.
32    pub order: KuramotoOrder,
33
34    /// Coupling coefficients `K_{i,j}` between the oscillators (nonnegative).
35    pub coupling_coeffs: DMatrix<f32>,
36
37    /// Damping coefficients `α_i` (nonnegative).
38    pub damping_coeffs: DVector<f32>,
39
40    /// Forcing parameter `P_i` (arbitrary sign).
41    ///
42    /// In the first-order case, these constant offsets are the *inherent
43    /// frequencies* of the oscillators, but in the second-order case, their
44    /// relationship to the inherent frequencies is more complicated; see
45    /// [Nishikawa & Motter
46    /// 2015](https://doi.org/10.1088/1367-2630/17/1/015012).
47    pub forcing_params: DVector<f32>,
48}
49
50impl KuramotoSystem {
51    /// Constructs a homogeneous Kuramoto system on the complete graph.
52    pub fn fully_connected_homogeneous(order: KuramotoOrder, forcing_params: DVector<f32>) -> Self {
53        let n = forcing_params.len();
54        let mut coupling_coeffs = DMatrix::from_element(n, n, 1.0);
55        for i in 0..n {
56            coupling_coeffs[(i, i)] = 0.0;
57        }
58        let damping_coeffs = DVector::from_element(n, 1.0);
59        Self {
60            order,
61            coupling_coeffs,
62            damping_coeffs,
63            forcing_params,
64        }
65    }
66}
67
68impl ODESystem for KuramotoSystem {
69    fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
70        match self.order {
71            KuramotoOrder::First => {
72                let n = x.len();
73                for i in 0..n {
74                    let mut rhs = self.forcing_params[i];
75                    for j in 0..n {
76                        rhs -= self.coupling_coeffs[(i, j)] * (x[i] - x[j]).sin();
77                    }
78                    dx[i] = rhs / self.damping_coeffs[i];
79                }
80            }
81            KuramotoOrder::Second => {
82                let n = x.len() / 2;
83                for i in 0..n {
84                    dx[i] = x[i + n];
85                    let mut rhs = self.forcing_params[i] - self.damping_coeffs[i] * x[i + n];
86                    for j in 0..n {
87                        rhs -= self.coupling_coeffs[(i, j)] * (x[i] - x[j]).sin();
88                    }
89                    dx[i + n] = rhs;
90                }
91            }
92        }
93    }
94}
95
96#[cfg(test)]
97mod tests {
98    use expect_test::expect;
99    use std::f32::consts::{FRAC_PI_2, FRAC_PI_4, TAU as TWO_PI};
100
101    use super::super::{ODEProblem, textplot_mapped_ode_result};
102    use super::*;
103
104    fn angular_distance(t1: f32, t2: f32) -> f32 {
105        let d1 = (t1 - t2) % TWO_PI;
106        let d2 = TWO_PI - d1;
107        d1.min(d2)
108    }
109
110    #[test]
111    fn first_order_kuramoto() {
112        let sys = KuramotoSystem::fully_connected_homogeneous(
113            KuramotoOrder::First,
114            DVector::from_element(3, 10.0 * TWO_PI),
115        );
116        let x0 = DVector::from_column_slice(&[-FRAC_PI_2, 0.0, FRAC_PI_4]);
117        let problem = ODEProblem::new(sys, x0).end_time(1.5);
118
119        // Check synchronization: convergence of angular distances to zero.
120        let result = problem.solve_rk4(0.001).unwrap();
121        let expected = expect![[r#"
122            ⡳⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 2.4
123            ⠄⠹⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
124            ⠂⠀⠘⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
125            ⡁⠀⠀⠘⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
126            ⠄⠀⠀⠀⠈⢧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
127            ⠂⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
128            ⣁⠀⠀⠀⠀⠀⠀⠳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
129            ⠌⠳⣄⠀⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
130            ⠂⠀⠘⢦⡀⠀⠀⠀⠀⠙⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
131            ⡁⠀⠀⠀⠙⢦⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
132            ⠄⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠳⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
133            ⠂⠀⠀⠀⠀⠀⠀⠈⠳⡄⠀⠀⠀⠙⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
134            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠘⢦⡀⠀⠀⠈⠳⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
135            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
136            ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦⣄⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
137            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⢤⡀⠀⠙⠲⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
138            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠲⢤⣀⠈⠙⠦⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
139            ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⠲⢤⣀⠉⠒⠦⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
140            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠓⠲⠤⢭⣙⣒⡦⠤⣄⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
141            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠙⠒⠚⠛⠻⠷⠶⠶⠶⠦⣤⣤⣤⣤⣤⡀
142            ⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠁ 0.0
143            0.0                                            1.5
144        "#]];
145        expected.assert_eq(&textplot_mapped_ode_result(
146            &problem,
147            &result,
148            |_i| true,
149            |x, i| angular_distance(x[i], x[0]),
150        ));
151    }
152
153    #[test]
154    fn second_order_kuramoto() {
155        let sys = KuramotoSystem::fully_connected_homogeneous(
156            KuramotoOrder::Second,
157            DVector::from_column_slice(&[1.0, -1.0, 1.0]),
158        );
159        let x0 = DVector::from_column_slice(&[-FRAC_PI_2, 0.0, FRAC_PI_4, 0.0, 0.0, 0.0]);
160        let problem = ODEProblem::new(sys, x0).end_time(10.0);
161
162        // These initial and parameter values are in the basin of attraction for
163        // the steady state with constant differences between phase angles.
164        let result = problem.solve_rk4(0.01).unwrap();
165        let expected = expect![[r#"
166            ⡙⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 2.4
167            ⠄⠈⢧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
168            ⠂⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
169            ⡁⠀⠀⢸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
170            ⠝⢦⠀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
171            ⠂⠈⣆⠀⢸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
172            ⡁⠀⠸⡀⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
173            ⠄⠀⠀⢧⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
174            ⠂⠀⠀⠸⡄⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
175            ⡁⠀⠀⠀⢇⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
176            ⠄⠀⠀⠀⢸⡀⠈⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
177            ⠂⠀⠀⠀⠀⡇⠀⢹⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣠⠞⠉⠀⠀⠈⠙⠲⢤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⣀⣀⣀⣀⣀⡀⠀⠀⠀
178            ⡉⠉⠉⠉⠉⢹⠉⠉⣏⠉⠉⠉⠉⠉⠉⠉⠉⢉⡟⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠙⠛⠻⠭⠭⠟⠛⠛⠋⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠁
179            ⠄⠀⠀⠀⠀⠈⡇⠀⠸⡄⠀⠀⠀⠀⠀⠀⣠⠏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
180            ⠂⠀⠀⠀⠀⠀⢱⠀⠀⠱⡀⠀⠀⠀⢀⡴⠃⠀⠀⠀⠀⣀⡤⠤⠤⢤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
181            ⡁⠀⠀⠀⠀⠀⠈⡆⠀⠀⠙⢤⣀⡤⠊⠀⠀⠀⠀⣠⠞⠁⠀⠀⠀⠀⠈⠙⠒⠤⣄⣀⠀⠀⠀⠀⠀⢀⣀⣀⣀⡤⠤⠤⠤⠤⠤⠤⠤⠤⠄
182            ⠄⠀⠀⠀⠀⠀⠀⢹⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡴⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠉⠉⠉⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
183            ⠂⠀⠀⠀⠀⠀⠀⠈⣇⠀⠀⠀⠀⠀⠀⠀⢀⠞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
184            ⡁⠀⠀⠀⠀⠀⠀⠀⠘⣆⠀⠀⠀⠀⠀⣠⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
185            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠘⢦⡀⠀⣠⠞⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
186            ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ -1.6
187            0.0                                           10.0
188        "#]];
189        expected.assert_eq(&textplot_mapped_ode_result(
190            &problem,
191            &result,
192            |i| i < 3,
193            |x, i| x[i] - x[0],
194        ));
195    }
196}
catlog/simulate/ode/kuramoto.rs

catlog/simulate/ode/
kuramoto.rs
