使用 MATLAB 的 ode45 和粗略的射击算法搜索 Mathieu 方程的周期解

计算科学 matlab 颂
2021-12-22 13:21:19

我正在使用 ODE45 对 Mathieu 方程进行数值模拟，并且我必须不断更改每个模拟的参数 delta 和 epsilon 以获得所需的周期性解。

以下是我正在运行的 MATLAB 代码：

function xdot=trial2(t,x) 
delta=0.1045;epsilon=0.0048685;
xdot=[x(2);(-delta-epsilon*cos(t))*x(1)-0.7*delta*abs(x(1))];

[t,x]=ode45('trial2',[0 10000000],[0;1]); 
hold on; plot(t,x(:,1),'r'); 
clear all;

如果我模拟 10000000 次，每次模拟需要 10 分钟，如果我将时间增加 5 倍，它会加倍。每次模拟后，我都会分析结果并相应地更改参数，以便更接近周期解。我必须为大约 30-40 个参数执行此操作，也许更多。我在具有 8gb 内存的 64 位桌面上运行它。该程序在 32 位桌面上运行内存不足。

1）有什么办法可以提高我的计算时间吗？我可能不得不从 10000000 到 50000000，甚至更多。2）你能帮我实现 parfor 循环，以便我可以同时运行更多参数的模拟。这也可能有一些帮助。3）如果我在比 i5-2400 CPU 更快的处理器上运行它会有帮助吗？

谢谢你。

2个回答

你的传播必须有多准确？您可以使用odeset函数来提高准确性（为 ode45 创建一组选项）。由于 ode45 是可变步长，因此精度对计算速度有很大影响。所需的精度越高，ode45 将采取的步骤越多。默认情况下，相对精度设置为 1e-3，绝对精度设置为 1e-6。

您还可以考虑使用其他语言（例如 C 或 Fortran）和Mex 函数或加载库（请参阅loadlibrary）进行集成。这肯定会导致巨大的加速，因为编译语言在这种类型的操作上比 Matlab 快得多。

编辑：添加到 Mex 功能的链接（对新用户的限制每个帖子只允许 2 个链接）

1）除非您在问题中提供的功能发生变化，否则功能的周期应该是恒定的。在这种情况下，我建议使用ode45定位“事件”的能力——例如过零（参见 文档odeset）。这样，您可以模拟更短的时间，并使用过零的准确估计来测量周期。这应该通过减少所需的求解器步骤数来减少您的计算时间。

2）您当然可以使用 parfor 循环来探索参数：

vfDelta = linspace(0, 1, 100);
vfEpsilon = linspace(0, 1, 100);
[mfD, mfE] = ndgrid(vfDelta, vfEpsilon);
mfParams = [mfD(:) mfE(:)];

parfor nParam = 1:size(mfParams, 1)
   % - Set up an anonymous function to pass to ode45, that contains your parameters
   fhSystem = @(t,x)trial2_params(t, x, mfParams(nParam, 1), mfParams(nParam, 2));

   % - Solve the ODE
   [t,x]=ode45(trial2, [0 1000], [0;1]);

   % - Do your analysis here...
   vfPeriod(nParam) = ...
end

function trial2_params(t, x, fDelta, fEpsilon)
   ...
end

3) 更快的 CPU 总是有帮助的……

其它你可能感兴趣的问题
上一篇土坡有限元法 下一篇我应该使用 bestman2 还是 dcache CLI（srmping、srmls 等）