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向量中？

我一直在编写我的代码，并使用events函数，我现在能够得到汽车的“齿轮”变化(作为事件)。现在，我遇到了一个与相同代码相关的新问题。想象一下，当我计算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将在每次t时计算加速度(y(1))和加速度的时间导数(y(2))。这似乎不正确.我认为初始值应该是初始位置(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) ) = [];

我要再说一遍:这是，而不是，我会推荐这个方法。只在测试或执行一些快速的n-脏东西时使用全局变量，但总是非常迅速地转向其他解决方案。而且，全局变量在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的每个评估按秒或更长的顺序进行时，这才是一个真正的问题，我怀疑这在您的设置中不是这样的。