# Instability of Kuramoto solutions with Runge Kutta 4th order method

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?

• Plots from top and down: 1) dt = 1e-1 2) dt = 1e-1 (so this is actually similar to plot 3), but on a shorter time scale) – Fabric Mar 17 '16 at 13:55
• In RK4, one sets $k_4= f(W+dt\cdot k_3, T + dt\cdot l_3)$ and $l_4= g(W+dt\cdot k_3, T + dt\cdot l_3)$, not $k_4= f(W+\frac12dt\cdot k_3, T + \frac12dt\cdot l_3)$ and $l_4= g(W+\frac12dt\cdot k_3, T + \frac12dt\cdot l_3)$ as written in your post. – Did Mar 17 '16 at 14:10
• I will edit that right away. I didn't mean to write that. I actually used the $k_4$ and $l_4$ that you wrote. But thanks! – Fabric Mar 17 '16 at 14:42
• Have you tried a different method, e.g. ode45? You may in fact need an implicit method to avoid instability. – Ian Mar 17 '16 at 14:48
• Haven't so far tried any other method. I was advised to use the RK4 method, but I vaguely remember something about RK4 and stiff equations not being a good match. What is ode45? – Fabric Mar 17 '16 at 14:58