3

MATLAB の微分方程式系で方程式の 1 つの位置を特定しようとしています。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

基本的に、ドキュメントを見て、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) が負の場合、ゼロ交差が 1 つしか表示されないことがあります。

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

    0.1085

ここに画像の説明を入力

于 2013-05-17T03:07:43.367 に答える
2

質問のタイトルが示すように、ODE の最大値を探している場合は、非常に近いです。微分方程式自体の根を使用して、これらの点を見つけます。つまり、導関数がゼロの場合です。これは、値がゼロ (またはその他の値) の解とは少し異なりますが、関連しています。問題は、指定value = Lorenz(t,x)していて、 のみに関心がある場合に ODE 関数がベクトルを返すことですx(3)。ただし、状態ベクトルにアクセスでき、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;

または、3 つの次元すべての最大値が必要な場合:

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

グローバルな最大値に関心がある場合は、出力を処理する必要がありますxm。または、システムが特定の振動動作をしている体制にいる場合は、 の最初の値に切り替えるか、単に見るisterminalだけでよい場合があります。1xm

最後に、無名関数を介してパラメーターを渡すことを検討してください。

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 に答える