I am currently trying to solve the Kuramoto model:
$\ddot{\theta_i} = P_i - \alpha\dot{\theta_i} + K \underset{i \neq j}{\sum}\sin(\theta_i - \theta_j)$
I split this second order differential equation into two ODEs:
$$\dot{\omega_i} = P_i - \alpha\omega_i + K \underset{i \neq j}{\sum}sin(\theta_i - \theta_j) := f(\omega_i, \theta_i)$$ $$\dot{\theta_i} = \omega_i := g(\omega_i, \theta_i)$$
Next I use the Runge Kutta 4th order method by approximating $f(\omega_i, \theta_i)$ and $g(\omega_i, \theta_i)$.
My concrete situation is a three node topology, that is $\omega_1, \omega_2, \omega_3$ and $\theta_1, \theta_2$ and $\theta_3$.
In this context, I consider P to be a 3-dim vector, W to be a 3-dim vector containing $\omega_i$ and T to b a 3-dim vector containing $\theta_i$.
My Runge kutta scheme then looks like:
$$k_1 = f(W, T)$$ $$l_1 = g(W,T)$$ $$k_2 = f(W+\frac{dt}{2}k_{1}, T + \frac{dt}{2}l_{1})$$ $$l_2 = g(W+\frac{dt}{2}k_{1}, T + \frac{dt}{2}l_{1})$$ $$k_3 = f(W+\frac{dt}{2}k_{2}, T + \frac{dt}{2}l_{2})$$ $$l_3 = g(W+\frac{dt}{2}k_{2}, T + \frac{dt}{2}l_{2})$$ $$k_4 = f(W+dtk_{3}, T + dtl_{3})$$ $$l_4 = g(W+dtk_{3}, T + dtl_{3})$$
Which gives me:
$$W(t+dt) = W(t) + \frac{dt}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$ $$T(t+dt) = T(t) + \frac{dt}{6}(l_1 + 2l_2 + 2l_3 + l_4)$$
Now when I run my program with a short pertubation $P(1) = -4.8$ from $t = 40$ to $t = 42$ with $dt = \frac{1}{10}$ I get a "expected" stability result. If I run my program with $P(1) = -4.8$ with $dt = \frac{1}{100}$ my solution completely blows up. Similar, if I run my program with $P(1) = -4.9$ and $dt = \frac{1}{10}$ my solution blows up aswell.
The plots are shown below. ON the Y-axis I plot $\theta_1 - \theta_2$ and on the X-axis I plot the time. Any idea what could cause this instability?
