1

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
%[Kuramoto Equation][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 zeor 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
%here I have only two oscillators from t = [0.0,0.5]

There is a small difference between the two matrices (which grows at larger time intervals). But unusually the total value calculated at each time (ie. each column) is the same. 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) and I think I am calling ode45 incorrectly, but I am not sure and haven't been able to figure out what is wrong for a few days. Any help would be appreciated.

James Kl
  • 177
  • 9
  • Perhaps this is germane: https://stackoverflow.com/questions/21859870/absolute-error-of-ode45-and-runge-kutta-methods-compared-with-analytical-solutio – duffymo Aug 28 '18 at 18:35
  • Not exactly as I am getting significantly different values between the two methods. – James Kl Aug 28 '18 at 21:07

2 Answers2

4

You should use the ode45 output. Runge Kutta as you have it implemented will eventually be unstable if you choose a step size that is too large. The entire point of ode45 is that it internally runs a Runge Kutta 4 and Runge Kutta 5 scheme. If the results of one integration step differ, then ode45 will reduce the time step until the results are comparable. Using the raw method like you are doing will obviously not do that.

Technically, things like ode45 are called "embedded Runge Kutta" methods. Here is one such method: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method They are efficient because the Runge Kutta methods of different order reuse a lot of the same function evaluations.

All this being said, you should find that if you reduce your time step enough, that the results are almost identical. The only reason that they differ is that ode45 is internally refining the time step when it detects that the solution may be inaccurate.

bremen_matt
  • 6,902
  • 7
  • 42
  • 90
  • Would that however explain why the difference consistently sums to 0? – James Kl Aug 28 '18 at 18:38
  • It is not clear to me from the question what point you are trying to make there. You are showing two rows of output. Are theses values at consecutive times? – bremen_matt Aug 28 '18 at 18:41
  • Also, if you suspect an error, the first thing I would do is to pull up a plot of the response using both methods. Are they substantively different? – bremen_matt Aug 28 '18 at 18:43
  • Sorry I do not know anything about your model. What do you mean when you say 2 different oscillators? Do you mean that you are running ode45 twice with different models? – bremen_matt Aug 28 '18 at 18:55
  • And what values does the Runge Kutta method return? – bremen_matt Aug 28 '18 at 19:01
  • I am not sure if that would help. The values don't look that different and both are reasonable. At time 200 the difference in the two methods is 60, which I think is too much (but I don't have enough context to tell). – James Kl Aug 28 '18 at 21:04
  • In your question, you give 2 different outputs of ode45. All I am saying is that you should add the output of Runge Kutta 4 for those same inputs, so that we can compare the outputs of both methods – bremen_matt Aug 29 '18 at 04:46
  • As to using the ode45 output: That would be the better result, only that in the presented code, the ode45 integration is for the wrong ODE, it integrates the ODE where its right side is always evaluated at the initial value, while the RK4 integration uses the ODE function consistently. – Lutz Lehmann Sep 08 '18 at 10:45
0

Did you actually use the line 

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

in your running code? Then the result here is definitely wrong. Use 

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

to get results that are at least related to the RK4 integration. That is, use the declared local variable/parameter of the anonymous function in the computation of its value. (Using u instead of y or theta to not reuse a variable name that is also used in the more global scope. Could use thetalocal instead if self-documenting variable names were desired.)

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • To clarify is u equal to theta(:,1)? – James Kl Sep 07 '18 at 04:08
  • No, why should it be? theta(:,1) is the initial value, the function you are passing has to be the ODE system function, which depends on the state passed in the parameters. – Lutz Lehmann Sep 07 '18 at 06:11
  • Sorry, could you clarify a little bit more. Specifically "use the declared local variable/parameter of the anonymous function in the computation of its value". Are you saying u should be symbolic? – James Kl Sep 07 '18 at 18:53
  • No, not symbolic. Read up on what anonymous functions resp. lambda expressions are. In your version, you are integrating a constant with ode45 while the RK4 integration integrates the ODE correctly. – Lutz Lehmann Sep 08 '18 at 10:43