我在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,它让我有一个巨大的速度,调用加速度计算(get acceleration):
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值抛出步进而不集成它?
谢谢你们。
最佳答案
首先:为什么微分方程的初始值是初始速度(speedInitial
)和初始加速度(accelerationInitial
)这意味着微分方程将计算加速度(car
)和加速度的时间导数jerk(y(1)
)这似乎不正确…我认为初始值应该是初始位置(y(2)
)和初始速度(t
)但是,我不认识你的模特,我可能错了。
现在,直接获取解决方案中的positionInitial
:不能,除非黑客入侵speedInitial
这也是合乎逻辑的;你也不能一直直接得到gear
,对吗这不是如何设置ode45
的。
我看到有两条出路:
全局变量
免责声明:不要使用这种方法这只是为了展示大多数人在第一次尝试时会做些什么。
您可以将dy
存储在全局变量中这可能是最少的编码量,但也是最不方便的结果:
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
的性质,这将得到一个时间数组gear
和相关的ode45
,其中包含中间点和/或被ts
拒绝的点所以,你必须在之后过滤:ts( ~ismember(ts, t) ) = [];
我再说一遍:这不是我推荐的方法只在测试或做一些快速n脏的事情时使用全局变量,但总是非常快速地转向其他解决方案此外,全局变量在每次(子)迭代
gear
时都会增长,这是不可接受的性能损失。最好使用下一种方法:
解决后呼叫
这对你的案子来说也不难,我建议你这样做首先,将微分方程修改如下,并按正常解:
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
这会让你得到所有的齿轮,汽车会有时间
ode45
。我在这里看到的唯一缺点是
ode45
的函数求值要比t
本身使用的函数求值多得多但这只是一个真正的问题,如果对car
的每次评估都需要几秒或更长的时间,我怀疑在您的设置中不是这样的。