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.