我不确定你是否能完全按照你的意愿去做,但是可以通过事件做很多事情。首先,这听起来像是某种数值计算首次通过时间 http://en.wikipedia.org/wiki/First-hitting-time_model(又名首次击中时间)。如果这些“粒子”是随机的,请停止并不要使用ode45
而是使用适合 SDE 的 pmethod。
据我所知,您可以拥有多少个事件函数 - 或者更确切地说,事件函数的维度(类似于 ODE 函数的维度) - 并且它们的数量与您拥有多少个 ODE 方程无关。 events 函数接收当前时间和当前状态向量。您可以使用其中的任何或所有元素来创建每个事件。您是正确的,更多的事件功能和更复杂的事件会减慢集成速度。性能还取决于检测事件的频率。如果你的每个粒子都到达“屋顶”,正如你所说的那样,并且只触发一个事件,那么这不会太糟糕。
在实现方面,这里是一个基于Matlab的简单例子ballode
例如,仅在垂直轴上模拟 N 个弹道粒子。有 N 个非终止事件来捕获每个粒子通过 y = 0 时的时间和速度。添加一个额外的终止事件来检查是否所有粒子都已通过 y = 0(如果我们知道这是哪个粒子)就像我们在这里所做的那样,我们可以将该事件设为终止事件)。
function eventsdemo
% Initial conditions for n balls
n = 10;
y0(2*n,1) = 0;
y0(n+1:end) = linspace(20,40,n);
% Specify events function
options = odeset('Events',@(t,y)efun(t,y,n));
% Integrate
[t,y,te,ye,ie] = ode45(@(t,y)f(t,y,n),[0 10],y0,options);
figure;
plot(t,y(:,1:n),'b',te(1:n),ye(sub2ind(size(ye),ie(1:n),(1:n).')),'r.');
function dydt = f(t,y,n)
% Differential equations for ballistic motion
dydt = [y(n+1:end);zeros(n,1)-9.8];
function [value,isterminal,direction] = efun(t,y,n)
% Last event checks that all balls have hit ground and terminates integration
yn = y(1:n);
value = [yn;all(yn < 0)];
zn = zeros(n,1);
isterminal = [zn;1];
direction = [zn-1;1];
在某些方面,这稍微低效,因为即使其中一些系统已经通过 y = 0,我们仍会继续模拟所有 N 个系统。但是,它很简单,并且输出数组是矩形的而不是参差不齐的。
我不清楚“耦合在一起”和需要粒子“停止”到底是什么意思。如果您需要做的不仅仅是记录事件数据,例如更改系统参数或以其他方式更改微分方程,那么您需要在每个事件后终止并重新启动积分。看着那(这ballode
示例(类型edit ballode
在 Matlab 命令窗口中)查看一些提高效率的建议。