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)。但是，我不知道你的型号，我可能是错的。

现在，直接进入gear解决方案：你不能，除非没有 hacking 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需要几秒钟或更长时间，我怀疑在您的设置中并非如此。