Modeling Kuramoto in Matlab

Question

I am trying to model Kuramoto ocillations in Matlab. $$ \frac{d\theta_i}{dt}=F_i(\theta)=\Omega_i+\frac{K}{N}\sum_{j=1}^N\sin(\theta_i-\theta_j),~~~i=1,...,N $$ I tried using ode45 to solve the system. I also saw someone else use the Runge-kutta method. I understand that ode45 uses the Runge-kutta method,however, the values I obtain from each are suspiciously different.

kuramoto= @(x,K,N,Omega)Omega+(K/N)*sum(sin(x*ones(1,N)-(ones(N,1)*x')))'
%Kuramoto is a model of N coupled ocilators (such as multiple radiowaves)
%The solution to the model is the phase of each ocilator


N = 2;
omega = rand(N,1);
theta(:,1) = 2*pi*randn(N,1);
t0 = theta(:,1);
[t,y] = ode45(@(t,y)kuramoto(theta(:,1),K,N,omega),tspan,t0);

%Runge-Kutta method
for j=1:iter
k1=kuramoto(theta(:,j),K,N,omega);
k2=kuramoto(theta(:,j)+0.5*h*k1,K,N,omega);
k3=kuramoto(theta(:,j)+0.5*h*k2,K,N,omega);           
k4=kuramoto(theta(:,j)+h*k3,K,N,omega);
theta(:, j+1)=theta(:,j)+(h/6)*(k1+2*k2+2*k3+k4);
end

Both methods output a matrix with N rows(where each row represents a different oscillator) and M columns (where M represents the solution at a given time) I have ode45 provide solutions form 0 to 0.5 at 0.1 intervals. To compare the methods I subtract the matrix obtained from Runge-Kutta with the matrix obtained using ode45. Ideally, the two should have the same values and the result should be a zero matrix but instead I get values such as:

0   -0.0003   -0.0012   -0.0027   -0.0048   -0.0076
0    0.0003    0.0012    0.0027    0.0048    0.0076

There is a small difference between the two matrices (which grows at larger time intervals). But unusually, the total theta calculated at each time (ie. each column) is the same between the two methods. This is consistent regardless of the number of oscillators.

I am unsure if this is a Math problem or programming problem (it's probably both). I am fairly confident in the implementation if the Rutta-kutta method as it is not my own code. Rather I believe I am missing something fundamental when I call ode45. I have been stuck for days and any help would be really appreciated.

Answer 1

One possible explanation, especially considering that the total theta is always in agreement between the two solvers, is that ode45 is a adaptive step solver, whereas the Runga-Kuta method (RK4) is a fixed step solver. Ode45 dynamically changes the time step it integrates over, based on the rate at which the solution is changing, in order to improve accuracy.

What this means it that the results of ode45 are spaced unevenly in time. By performing an immediate elementwise comparison between the results of ode45 and RK4, you are comparing solution values at between different times. In order to make a useful comparison, you either need to pass the desired output time vector to ode45 when you call it (this does not affect the calculation process, only the output spacing), or to interpolate the results of ode45 at the desired output times.

I apologize if you already knew this. It is unclear in your question whether you have accounted for the uneven timing, and I don't yet have the reputation to comment for clarification. I know I was caught by this behavior when first working with ODEs in Matlab.

Answer 2

The observed error results from using the wrong ODE. What you want to solve is $\dot\theta=F(\theta)$. But by calling the solver slightly wrong with

[t,y] = ode45(@(t,y)kuramoto(theta(:,1),K,N,omega),tspan,t0);

you actually solve $\dot\theta=F(\theta(0))$. The expression in the place of the function argument,

@(t,y)kuramoto(theta(:,1),K,N,omega)

defines an anonymous function or lambda expression. The local arguments of that function are given in @(t,y), all the remaining variables are passed as constants. The result is that this ODE function neither depends on time nor on the state y, it returns a constant slope, giving linear functions as solutions.

However, the ODE function $F$ in the original task depends on the state, y in @(t,y) needs to be passed to the function computing the derivatives from it, $(t,y)\mapsto F(y)$. Also, it may not matter, but re-using the same letter in two different roles in the same line is at least confusing. With

Omega = rand(N,1);
theta0(:,1) = 2*pi*randn(N,1);
[t,y] = ode45(@(t,u)kuramoto(u,K,N,Omega),tspan,theta0);

you should get correct results are close to the RK4 solution.

---

PS: That the difference of the solutions adds up to zero is due to the fact that the derivative vectors add up to the constant sum(Omega) in both cases, so that the time evolution of the sums over the states are the same linear function, thus giving zero in their difference. If the state vector only has 2 elements, the elements of the difference vector thus have to be opposite, as observed.