4

オンラインでこの投稿を読んでいたところ、「if ステートメント」と「abs()」関数を使用すると、MATLAB の可変ステップ ODE ソルバー (ODE45 など) に悪影響を与える可能性があると述べられていました。OP によると、時間ステップに大きな影響を与え (必要な時間ステップが低すぎる)、微分方程式が最終的に積分されたときに悪い結果が生じる可能性があります。これが本当かどうか、もしそうなら、なぜだろうと思っていました。また、修正ステップ ソルバーに頼らずに、この問題を軽減するにはどうすればよいでしょうか。私が何を意味するかについて、以下のコード例を示しました。

function [Z,Y] = sauters(We,Re,rhos,nu_G,Uinj,Dinj,theta,ts,SMDs0,Uzs0,...
Uts0,Vzs0,zspan,K)

Y0 = [SMDs0;Uzs0;Uts0;Vzs0]; %Initial Conditions
options = odeset('RelTol',1e-7,'AbsTol',1e-7); %Tolerance Levels
[Z,Y] = ode45(@func,zspan,Y0,options);

function DY = func(z,y)

    DY = zeros(4,1);

    %Calculate Local Droplet Reynolds Numbers
    Rez = y(1)*abs(y(2)-y(4))*Dinj*Uinj/nu_G;
    Ret = y(1)*abs(y(3))*Dinj*Uinj/nu_G;

    %Calculate Droplet Drag Coefficient
    Cdz = dragcof(Rez);
    Cdt = dragcof(Ret);

    %Calculate Total Relative Velocity and Droplet Reynolds Number
    Utot = sqrt((y(2)-y(4))^2 + y(3)^2);
    Red = y(1)*abs(Utot)*Dinj*Uinj/nu_G;

    %Calculate Derivatives
    %SMD
    if(Red > 1)
        DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...
            Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...
            (We/Re)*K*(Red^0.5)*Utot*Utot/y(2);
    elseif(Red < 1)
        DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...
        Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...
        (We/Re)*K*(Red)*Utot*Utot/y(2);
    end
    %Axial Droplet Velocity
    DY(2) = -(3/4)*rhos*(Cdz/y(1))*Utot*(1 - y(4)/y(2));
    %Tangential Droplet Velocity
    DY(3) = -(3/4)*rhos*(Cdt/y(1))*Utot*(y(3)/y(2));
    %Axial Gas Velocity
    DY(4) = (3/8)*((ts - ts^2)/(z^2))*(cos(theta)/(tan(theta)^2))*...
        (Cdz/y(1))*(Utot/y(4))*(1 - y(4)/y(2)) - y(4)/z;

end

end

関数「dragcof」は次のように指定されます。

function Cd = dragcof(Re)    
if(Re <= 0.01)

    Cd = (0.1875) + (24.0/Re);

elseif(Re > 0.01 && Re <= 260.0)

    Cd = (24.0/Re)*(1.0 + 0.1315*Re^(0.32 - 0.05*log10(Re)));

else

    Cd = (24.0/Re)*(1.0 + 0.1935*Re^0.6305);
end
end
4

2 に答える 2

4

これは、ifステートメント、モジュラス演算 ( abs())、または単位ステップ関数、ディラック デルタなどを使用して計算される導関数が、解またはその導関数の値に不連続性を導入し、ねじれやジャンプが発生するためです。 、変曲点など

これは、ODE の解が、関連する時点で動作が完全に変化することを意味します。可変ステップ積分器が行うことは

  • これを検出
  • 「問題点」を超えて直接情報を使用できないことを認識する
  • ステップを減らし、問題点が精度の要求を満たすまで、上から繰り返します

したがって、多くの失敗したステップが発生し、問題点の近くでステップ サイズが縮小され、全体的な統合時間に悪影響を及ぼします。

ただし、可変ステップ積分器は引き続き良好な結果を生成します。一定ステップ積分器は、そもそもそのような問題を検出できないため (エラー推定がないため)、この種の問題の良い解決策ではありません。

あなたができることは、問題を複数の部分に単純に分割することです。変更が発生する時点が事前にわかっている場合は、前の積分の出力を次の積分の初期値として使用するたびに、間隔ごとに新しい積分を開始するだけです。

問題がどこにあるのか事前にわからない場合は、イベント関数と呼ばれる Matlab の ODE ソルバーの非常に優れた機能を使用できます(ドキュメントを参照してください)。Matlab のソルバーの 1 つにイベント (導関数の符号の変化、ifステートメントの条件の変化など) を検出させ、そのようなイベントが検出されたときに統合を終了します。次に、前回と同様に前回の積分の初期条件を使用して、最後の時間に達するまで、新しい積分を開始します。

Matlab はイベントの位置を正確に検出しようとするため、この方法でも全体の実行時間に若干のペナルティが発生します。ただし、実行時間と結果の精度の両方に関しては、やみくもに統合を実行するよりもはるかに優れています。

于 2013-05-16T09:16:19.577 に答える