3

我试图在 MATLAB 的微分方程系统中找到一个方程的位置。我正在尝试使用 odeset 的事件属性。如何在我的函数中挑选出特定的方程？

options = odeset('Events',@event);
[t x tm xm ie] = ode45(@Lorenz,[0 200],I,options);


function X = Lorenz(t,x)
r = 15;
sigma = 10;
b = 8/3;
X(1,1) = sigma*(x(2,1)-x(1,1));
X(2,1) = r*(x(1,1)) - x(2,1) -x(1,1)*x(3,1);
X(3,1) = x(1,1)*x(2,1) - b*x(3,1);
end

function [value,isterminal,direction] = event(t,x)
value  = Lorenz(t,x); %I can't specify X(3,1) here
isterminal = 0;
direction = -1;
end

特别是我试图在 X(3,1) = 0 时记录。

谢谢

4

2 回答 2

2

如果您正在寻找 ODE 的最大值，正如您的问题标题所示，那么您非常接近。您正在使用微分方程本身的根来找到这些点，即当导数为零时。这与具有零（或其他）值但相关的解决方案略有不同。问题是您正在指定value = Lorenz(t,x)并且 ODE 函数在您只对x(3). 但是您可以访问状态向量并拥有三种选择。

最简单的：

function [value,isterminal,direction] = event(t,x)
b = 8/3;
value = x(1)*x(2)-b*x(3); % The third equation from Lorenz(t,x) 
isterminal = 0;
direction = -1;

或者，效率较低：

function [value,isterminal,direction] = event(t,x)
y = Lorenz(t,x); % Evaluate all three equations to get third one
value = y(3);
isterminal = 0;
direction = -1;

或者，如果您想要所有三个维度的最大值：

function [value,isterminal,direction] = event(t,x)
value = Lorenz(t,x);
isterminal = [0;0;0];
direction = [-1;-1;-1];

如果您对全局最大值感兴趣，那么您需要处理输出，xm. 或者，如果您处于系统具有某些振荡行为的状态，那么您可以切换isterminal到1或只查看xm.

最后，您可以考虑通过匿名函数传递参数：

r = 15;
sigma = 10;
b = 8/3;
f = @(t,x)Lorenz(t,x,r,sigma,b);

I = [1;5;10];
options = odeset('Events',@(t,x)event(t,x,b));

[t,x,tm,xm,ie] = ode45(f,[0;10],I,options);

和：

function X = Lorenz(t,x,r,sigma,b)
...

function [value,isterminal,direction] = event(t,x,b)
...
于 2013-05-19T20:10:57.723 回答
2

基本上，查看文档，如果您有兴趣查看 x(3) = 0 的时间，那么您需要重写您的事件函数：

function [value,isterminal,direction] = event(t,x)
value  = x(3); %I can't specify X(3,1) here --> why not?? Lorenz(t,x) is going to return the differential. That's not what you want
isterminal = 0;
direction = 0; %The desired directionality should be "Detect all zero crossings"
end

现在我不知道你是如何定义I的

[t x tm xm ie] = ode45(@Lorenz,[0 200],I,options);

但是您对方程的解在几个点附近非常稳定，如果 x(3) 在零时为负，您可能只会看到一个零交叉。

[t x tm xm ie] = ode45(@khal,[0 5],[10 -1 -20],options);
tm =

    0.1085

在此处输入图像描述

于 2013-05-17T03:07:43.367 回答