0

MATLAB で ode45 を使用して一連の ODE を実行していますが、後で使用するために変数の 1 つ (導関数ではない) を保存する必要があります。関数 'assignin' を使用して、ベース ワークスペースに一時変数を割り当て、各ステップで更新しています。これは機能しているように見えますが、配列のサイズが ode45 から取得した解ベクトルのサイズと一致しません。たとえば、次のネストされた関数があります。

function [Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0)

options = odeset('RelTol',1e-7,'AbsTol',1e-7);
[Z,Y] = ode45(@momentum,zspan,Y0,options);

function DY = momentum(z,y)

    DY = zeros(4,1);

    %Entrained Total Velocity
    Ve = sin(theta)*(y(4));

    %Total Relative Velocity
    Urs = sqrt((y(1) - y(4))^2 + (y(2) - Ve*cos(theta))^2 + (y(3))^2);

    %Coefficients
    PSI = K*Urs/y(1);
    PHI = P*Urs/y(1);

    %Liquid Axial Velocity
    DY(1) = PSI*sign(y(1) - y(4))*(1 + (1/6)*(abs(y(1) - y(4))*G)^(2/3));

    %Liquid Radial Velocity
    DY(2) = PSI*sign(y(2) - Ve*cos(theta))*(1 + (1/6)*(abs(y(2) - ...
        Ve*cos(theta))*G)^(2/3));

    %Liquid Tangential Velocity
    DY(3) = PSI*sign(y(3))*(1 + (1/6)*(abs(y(3))*G)^(2/3));

    %Gaseous Axial Velocity
    DY(4) = (1/z/y(4))*((PHI/z)*sign(y(1) - y(4))*(1 + ...
        (1/6)*(abs(y(1) - y(4))*G)^(2/3)) + Ve*Ve - y(4)*y(4));

    assignin('base','Ve_step',Ve);
    evalin('base','Ve_out(end+1) = Ve_step');
end

end

上記のコードでは、theta (ラジアン)、K (負の値)、P、および G は定数であり、この例では、任意の値を使用できます。Zspan は ODE ソルバーの積分時間ステップであり、Y0 は初期条件ベクトル (4x1) です。繰り返しますが、この例のために、これらは任意の妥当な値を取ることができます。メインファイルでは、関数は次のように呼び出されます。

Ve_out = 0;
[Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0);
Ve_out = Ve_out(2:end);

この方法は MATLAB から問題なく動作しますが、問題は Ve_out のサイズが Z または Y のサイズと同じではないことです。この理由は、MATLAB がそのアルゴリズムに対して ODE 関数を複数回呼び出すためです。したがって、解決策は次のとおりです。 Ve_out よりわずかに小さくなります。am304 が示唆したように、ODE 関数に DY = モメンタム(Z,Y) などの Z および Y ベクトルを与えることで DY を単純に計算できますが、これを「assignin」(または同様の方法) で機能させる必要があります。この問題のバージョンでは、DY と Ve の間に暗黙の依存関係があり、反復ごとに DY を計算するには計算コストがかかりすぎます (この問題を何度も反復して実行します)。

4

2 に答える 2

3

それでは、SSCCE の簡単な例から始めましょう。

function [Z,Y] = khan

options = odeset('RelTol',1e-7,'AbsTol',1e-7);
[Z,Y] = ode45(@momentum,[0 12],[0 0],options);
end

function Dy = momentum(z,y)
Dy = [0 0]';

Dy(1) = 3*y(1) + 2* y(2) - 2;
Dy(2) = y(1) - y(2);

Ve = Dy(1)+ y(2);

    assignin('base','Ve_step',Ve);
    evalin('base','Ve_out(end+1) = Ve_step;');

    assignin('base','T_step',z);
    evalin('base','T_out(end+1) = T_step;');
end

コマンドラインとして実行[Z,Y] = khanすることで、問題を示す完全な機能コードを得ることができ、頭痛の種は一切ありません。これに対する私の忍耐力は尽きました。生きて学びましょう。

これは機能しているように見えますが、配列のサイズが ode45 から取得した解ベクトルのサイズと一致しません。

時間変数を抽出するコードに 2 行追加したことに注意してください。コマンド プロンプトから次のコマンドを実行するだけで、何が起こっているのかを理解できます。

Ve_out = [];
T_out = [];
[Z,Y] = khan;
size (Z)
size (T_out)
size (Ve_out)
plot (diff(T_out))

ans =
   109     1

ans =   
     1   163

ans =
     1   163

ここに画像の説明を入力

基本的に、ode45 は反復アルゴリズムです。つまり、定期的にコースが正しいことを意味します (これが、定期的に diff(T) = 0 が表示される理由です)。アルゴリズムに自分のやりたいことを強制することはできません。それを受け入れなければなりません。

したがって、選択肢は次のとおり
です。 1.固定ステップ アルゴリズムを使用します
。 2. ode45 アルゴリズムが汚れた作業を行った後に、必要なものを再現する関数呼び出しを行います。(am304 の解決策)
3. 時間関数を使用してデータを収集し、すべてのアルゴリズムを解析して余分なデータを削除します。

于 2013-05-22T18:47:47.913 に答える