使用ODE45查找函數的最大值

Question

我試圖在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。

謝謝

Answer 1

基本上，看着documentation，如果你有興趣看當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 

現在我不知道你是怎麼在

定義

[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 

Answer 2

如果你正在尋找一個ODE的最大值，因爲你的問題的標題所示，那麼你是非常接近的。您正在使用微分方程本身的根來找到這些點，即導數爲零時。這與具有零（或其他）值的解決方案略有不同，但相關。問題在於你指定value = Lorenz(t,x)，而當你只對x(3)感興趣時，ODE函數返回一個向量。但是你可以訪問狀態向量並有三種選擇。

最簡單的：

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]; 

如果你有興趣全球最大值，那麼您將需要處理輸出，​​。或者，如果您處於系統具有某些振盪行爲的制度，那麼您可能只需將isterminal切換爲1，或者只需查看​​中的第一個值即可。

最後，您可以考慮通過匿名函數傳遞你的參數：

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) 
...