3

我正在 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 值抛出步进而不集成它？

多谢你们。

4

1 回答 1

5

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

现在，直接进入gear解决方案：你不能，不是没有黑客攻击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，这是不可接受的性能损失。

最好使用下一种方法：

解决后调用

这对你的情况也不是太难，我建议你走的路。首先，修改微分方程如下，并正常求解：

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。

car我在这里可以看到的唯一缺点是，您将拥有比单独使用更多的功能评估ode45。但这只是一个真正的问题，如果每次评估都car需要几秒钟或更长时间，我怀疑您的设置不是这种情况。

于 2012-11-27T08:13:43.137 回答