matlab ode45 检索参数

Question

我正在 Matlab 中试验 ode45。我已经学会了如何将参数传递给 ode 函数，但我仍然有一个问题。假设我想计算汽车的轨迹(速度曲线)并且我有一个函数，例如getAcceleration，它为我提供了汽车的加速度以及正确的档位:[acceleration, gear] = getAcceleration(speed,modelStructure) 其中 modelStructure 代表汽车的型号。

ode 函数将是:

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);

然后我这样调用Ode45积分器:

tInit = 0;
tEnd  = 5,
[t,y] = ode45(@car,[tInit tEnd], [speedInitial,accelerationInitial],options,modelStructure);

问题是:如何获取存储齿轮的向量？我应该有类似 [t,y,gear]=ode45(....) 的东西还是 gear 应该在 y 向量中？

---

我一直在编写我的代码并使用事件函数，我现在能够获得汽车“齿轮”的变化(作为事件)。 现在我有一个与相同代码相关的新问题。 想象一下，当我评估 de 'dy' 向量时，我能够得到一个进一步的值 Z，这让我能够大幅加快调用加速度计算 (getAcceleration) 的速度:

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),Z(t)] = getAcceleration(y(1),modelStructure,Z(t-1)); 

并假设我也能够在初始条件下计算 Z。问题是我无法计算 Z 导数。

有没有办法在不集成的情况下传递 Z 值抛出步进？

谢谢你们。

Answer 1

首先:为什么微分方程的初始值是初始速度 (speedInitial) 和初始加速度 (accelerationInitial)？这意味着微分方程 car 将计算加速度 (y(1)) 和加加速度 (y(2))，加速度的时间导数，每次 t。这似乎不正确...我会说初始值应该是初始位置 (positionInitial) 和初始速度 (speedInitial)。但是，我不知道你的模型，我可能是错的。

现在，直接在解决方案中获取装备:你不能，除非破解ode45。这也是合乎逻辑的；你也不能一直直接得到 dy，对吧？ ode45 不是这样设置的。

我在这里看到了两条出路:

全局变量

免责声明:请勿使用此方法。这只是为了展示大多数人在第一次尝试时会做什么。

您可以将gear 存储在全局变量中。这可能是最少的编码，但也是最不方便的结果:

global ts gear ii

ii    = 1;
tInit = 0;
tEnd  = 5,
[t,y] = ode45(...
    @(t,y) car(t,y,modelStructure), ...
    [tInit tEnd], ...
    [speedInitial, accelerationInitial], options);

...

function [dy] = car(t,y,modelStructure)
global ts gear ii

dy    = zeros(2,1);
dy(1) = y(2);
[dy(2),gear(ii)] = getAcceleration(y(1),modelStructure);

ts(ii) = t;
ii = ii + 1;

但是，由于 ode45 的性质，这将为您提供一组时间 ts 和相关的 gear，其中包含中间点和/或被 ode45 拒绝的点。所以，你必须在之后过滤那些:

ts( ~ismember(ts, t) ) = [];

我再说一遍:这不是我推荐的方法。只在测试或做一些快速而肮脏的事情时使用全局变量，但总是很快转向其他解决方案。此外，全局变量在 ode45 的每次(子)迭代中增长，这是 Not Acceptable 性能损失。

最好使用下一个方法:

解决后调用

这对您的情况来说也不难，我建议您采用这种方式。首先，修改微分方程如下，正常求解:

tInit = 0;
tEnd  = 5,
[t,y] = ode45(...
    @(t,y) car(t,y,modelStructure), ...
    [tInit tEnd], ...
    [speedInitial, accelerationInitial], options);

...

function [dy, gear] = car(t,y,modelStructure)    

dy    = [0;0];
dy(1) = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);

然后在 ode45 完成后，执行此操作:

gear = zeros(size(t));
for ii = 1:numel(t)
    [~, gear(ii)] = car(t(ii), y(ii,:).', modelStructure); 
end

这将使您获得汽车在 t 时可能拥有的所有档位。

我在这里看到的唯一缺点是，与 ode45 本身相比，您对 car 的函数求值要多得多。但这只是一个真正的问题，如果 car 的每次评估都需要几秒或更长时间，我怀疑您的设置不是这种情况。