我正在尝试在Matlab中模拟仓本的OCILLATIONS。我尝试使用 ode45 来解决系统问题。我还看到其他人使用Runge-kutta方法。我知道 ode45 使用 Runge-kutta 方法，但是，我从每个方法获得的值都可疑地不同。
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
两种方法都输出一个矩阵，其中包含 N 行(其中每行代表不同的振荡器(和 M 列(其中 M 代表给定时间的解(，我让 ode45 以 0.1 的间隔提供从 0 到 0.5 的解决方案。为了比较这些方法，我将从Runge-Kutta获得的矩阵与使用ode45获得的矩阵减去。理想情况下，两者应该具有相同的值，结果应该是一个 zeor 矩阵，但我得到
的值如下：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]
两个矩阵之间存在很小的差异(以较大的时间间隔增长(。但不同寻常的是，每次(即每列(计算的总值是相同的。无论振荡器的数量如何，这都是一致的。
我不确定这是数学问题还是编程问题(可能两者兼而有之(，我想我错误地调用了 ode45，但我不确定并且几天来一直无法弄清楚出了什么问题。任何帮助将不胜感激。
您应该使用 ode45 输出。如果您选择的步长太大，则实现的 runge Kutta 最终将不稳定。ode45的全部意义在于它在内部运行Runge Kutta 4和Runge Kutta 5方案。如果一个积分步骤的结果不同，则 ode45 将减少时间步长，直到结果具有可比性。像您正在做的那样使用原始方法显然不会这样做。
从技术上讲，像ode45这样的东西被称为"嵌入式Runge Kutta"方法。这是一种这样的方法：https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method 它们之所以高效，是因为不同阶次的 Runge Kutta 方法重用了许多相同的函数计算。
综上所述，您应该发现，如果您足够减少时间步长，结果几乎相同。它们不同的唯一原因是 ode45 在检测到解决方案可能不准确时在内部优化时间步长。
你真的用过这条线吗
[t,y] = ode45(@(t,y)kuramoto(theta(:,1),K,N,omega),tspan,t0);
在您的运行代码中？那么这里的结果肯定是错误的。用
[t,y] = ode45(@(t,u)kuramoto(u,K,N,omega),tspan,t0);
以获得至少与 RK4 集成相关的结果。也就是说，在计算其值时使用匿名函数的声明局部变量/参数。(使用u而不是y或theta不重用也在更全局范围内使用的变量名。如果需要自记录变量名称，可以使用thetalocal。
PS：差和为零是由于导数向量的总和为零，因此状态向量上的和是一个常数，而不管在应用方法时犯了什么错误。因此，如果您从自身中减去相同的常数，则最终再次为零。如果状态向量只有 2 个元素，则差分向量的元素必须相反。