0

私の教科書では、ODEのシステムのイベントの場所を指定するときに使用することになっている関数の例に遭遇しました。関数の例は次のとおりです。

function [value, isterminal, dircn] = proj(t,z,flag);
g = 9.81;
if nargin < 3 | isempty(flag)
    value = [z(2); 0; z(4); -g];
else
    switch flag
        case 'events'
            value = z(3);
            isterminal = 1;
            dircn = -1;
        otherwise
            error('function not programmed for this event');
    end
end

ここに私が理解していない論理の一部があります。「イベント」オプションをアクティブにしてから、ode45を実行するとします。では、ode45は実際に連立方程式（上記の関数でとして指定されているvalue = [z(2); 0; z(4); -g];）をどのように読み取ることができますか？もちろん、tspanと初期条件を指定した後、上記の関数に基づいてode45を実行しましたが、これは魅力のように機能します。しかし、上記のスクリプトの「if」部分にのみ表示されている場合に、ode45がシステムを正しく読み取る方法がわかりません。

誰かがここで論理を説明することができれば、私はそれを大いに感謝します！

4

3 に答える 3

1

答えはifコマンドにあります：

if nargin < 3 | isempty(flag)
    value = [z(2); 0; z(4); -g];
else

引数の数が3未満の場合、または変数flagが空の場合は、変数valueをに設定します[z(2); 0; z(4); -g]。それ以外の場合、変数flagが、の場合は変数を'events'に設定し、そうでない場合はエラーを報告します。したがって、この関数は常に変数またはに何らかの戻り値を割り当て、コマンドを使用してエラーを報告します。valuez(3)flag'events'valueerror

于 2012-05-15T22:40:45.330 に答える
1

さて、私はいくつかの部分を説明できると思います。上にも書きましたが、次元がvalue変わるのは不思議です。

状態空間と変数の名前を考えると、2 次元の動きのように見えます。

フラグがない場合、状態空間は次のようになります。

  1. 水平位置 (x)
  2. 水平速度 (vx)
  3. 垂直位置 (y)
  4. 垂直速度 (vy)

訂正 ode は'events'指定すると送信できるようです。したがって、関数は状態空間の 3 番目のコンポーネントを出力します。それを説明しているこのサイトを見てください：http://www.mathworks.de/help/techdoc/math/f1-662913.html

事前に指定しないと、ode'events'は関数に送信できないため、この部分は呼び出されません。

ただし、関数はとにかく機能しません。導関数は状態空間 (4x1) と同じ次元である必要があるためです。しかし、1x1しかありません。

しかし、「イベントの場所を指定する」とは何を意味するのかよくわかりません。そこに秘密が隠されているのかもしれません。

創造性があれば、関数を使用して状態空間の 3 番目のコンポーネントを抽出できると思います。

于 2012-05-16T08:49:37.943 に答える
0

上記の関数を記述した後の処理方法について、もう少し情報を補足できます。次に、次のように定義します。

tspan = [0 6];
z0 = [0, 5*cos(pi/4), 0, 5*sin(pi/4)];
options = odeset('events','on');
[t y] = ode42('proj',tspan,z0,options)

これを記述すると、出力は 5 列になりtspanます。ゼロに達すると、「イベント」が発生するため、計算は終了します。z(1)z(2)z(3)z(4)z(3)

一方、このoptions = odeset('events','on')行を含めず、単純に次のように記述した場合:

[t y] = ode42('proj',tspan,z0)

計算は範囲全体に対して実行されtspanます。

それでも、ode42「イベント」をアクティブにしたときにすべての出力ベクトルを計算する方法をまだ確認できていません。MatLab が「proj」関数で「else」ステートメントのみを実行する必要があることが論理的に見えるからです。関数のこの部分には、微分方程式の実際のシステムは含まれていません。

于 2012-05-16T11:05:44.443 に答える