matlab - Matlab:ode45 和 4 阶 Runge-Kutta 方法产生不同的值
问题描述
我正在尝试在 Matlab 中模拟仓本振荡。我尝试使用 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 这样的东西被称为“嵌入式龙格库塔”方法。这是一种这样的方法:https ://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method 它们很有效,因为不同顺序的 Runge Kutta 方法重用了许多相同的函数评估.
说了这么多,您应该会发现,如果您将时间步长减少到足够多,那么结果几乎是相同的。它们不同的唯一原因是 ode45 在检测到解决方案可能不准确时会在内部优化时间步长。
推荐阅读
- java - 有什么方法可以读取在 GRID 机器(Selenium-Java)中下载的 csv 文件
- azure-eventhub - 事件中心捕获数据积压
- java - linux环境下无法解压zip文件
- javascript - 只选择一个复选框来切换一组 div
- ionic4 - 如何在侧面菜单离子 4 中获取数据
- excel - 如何仅从特定打开 Excel 文件并从每个 Excel 文件中复制特定单元格值
- vba - 制作打印的 Word 文档的副本并将其保存在文件夹中
- python - Python - 类装饰器 - 运行时继承
- php - TCPDF:如果打印 pdf,如何更新数据库打印状态
- typescript - 如果不能动态更改对象的类怎么办?