# Modeling Kuramoto in Matlab

I am trying to model Kuramoto ocillations in Matlab. 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.

• Ode45 is not RK4, it is Dormand-Price which is RK-type but not the same method. – Ian Aug 28 '18 at 22:59
• Just noticed that in the documentation for ode45. But would that significantly change the solution each method gets? – James Kl Aug 29 '18 at 3:40
• What is tspan? Are you sure you forced ode45 to give evenly spaced output? – Ian Aug 29 '18 at 4:09

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.

• I set the time interval which ode45 integrates over to be the same as the interval which RK4 integrates over. I was able to confirm it too by looking at the value of t. – James Kl Aug 29 '18 at 3:39

The observed error results from using the wrong ODE by calling the solver slightly wrong with

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


The state y in @(t,y) needs to be passed to the function computing the derivatives from it. Also, it may or may not matter, but re-using the same letter in two different roles in the same line is at least confusing. With

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


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

PS: That the difference sums to zero is due to the fact that the derivative vectors sum to zero, so that the sum over the state vector is a constant regardless of the errors committed in applying the methods. So if you subtract the same constant from itself you end up again with zero. If the state vector only has 2 elements, the elements of the difference vector thus have to be opposite.