非負の解を与えるように制約された遅延微分方程式（DDE）システムを解く

Question

ode45は、解を非負に制限するNonNegativeというパラメータを持っています。 They even wrote a paper about how this method worksと、負の値になるたびにy_iを0に設定するのと同じように愚かではないことが一般的にはうまくいきません。

現在、MATLABには遅延微分方程式を解くためのdde23がありますが、この積分器のパラメータはNonNegativeです。

残念ながら、NonNegativeを有効にしてode45を使用して解決される既存のODEシステムに遅延を追加することを任されます。

どのように進めるべきですか？

EDIT：

私はこれが役に立つかどうかわからないんだけど、...

私のシステムのDDE部分は、基本的には次のようになります。

dx = 1/(1+y*z) - x; 
dy = (y*z)^2/(1+(y*z)^2) - y; 
dz = X - z; 

X（資本金第3の式における文字変数）は、xの遅延バージョンである。次に、私はxとzの方程式に2つの項を加えて、組み合わせたシステムをすべて統合することで、このDDEシステムを既存の（そしてより大きな）ODEシステムにリンクしています。

Answer 1

難しい問題があり、ワンステップ解決策があるかどうかはわかりません。私は、代わりの答えを提供しようと思っている人には名誉を提供することをもっとうれしく思います。

遅延の長さによっては、方程式を複数回実行し、各反復で古い値のxを最新の更新に渡すことがあります。

たとえば、遅延が1時間であるとします。最初の1時間で、NonNegativeがフラグされたode45を実行します。 TimeパラメータとともにNew Matrixに値を格納し、アルゴリズムを再度実行します。今、洗いすすぎ、そして繰り返し古いソリューションマトリックスと昔行列

dx = 1/(1+y*z) - x; 
dy = (y*z)^2/(1+(y*z)^2) - y; 
tindex = find(told>t,1) -1 % find the upper index which best approximates t 
X = xold(tindex) + (xold(tindex+1)-xold(tindex))*(t-told(tindex))/(told(tindex+1)-told(tindex)) % or interpolation method of your choosing 
dz = X - z; 

：この時間は、あなたが2つの入力パラメータを追加してください。ここで、Xは、例3のように、準時間依存の項であるode45であることに注意してください。

Answer 2

私は最近、私のいくつかのコードでこの問題を抱えていました。「最も簡単なのソリューションは、ソリューションが0に達すると、まず、あなたは は、このように（x == 0ケースが評価された予告）

if x > 0 
    dx = 1/(1+y*z) - x; 
else % if x <= 0 
    dx = 0; 
end 

に

dx = 1/(1+y*z) - x; 

を変更する、0でそれを維持する必要があり、次のことを行うことです

dxTmp = 1/(1+y*z) - x; 
if x > 0 
    dx = dxTmp; 
elseif dxTmp > 0 
    % x becomes positive again 
    dx = dxTmp; 
else 
    dx = 0; 
end 

注ただし、これは一次導関数にジャンプ不連続部を作成すること（それが0になることはありませんできる理由に応じて）、WHIに

それともDDEソルバーがこの点の近くにt - delayを持つポイントに達すると、この不連続点が正確な場所を知らない限り、それを非常に効率的には解決しません（通常はMatlabにどこにあるかを伝える余分なオプションを使用しますが、以下の手順に従うと、それは必要ありません）。

この不連続箇所を特定するには、というDDEオプションを使用する必要があります（「イベントの場所プロパティ」までスクロールしてください）。examplesを見ることもできます。これらの例の1つは、値はODEでは使用できません - ODEとDDEのイベントはほとんど同じです）。基本的に、イベントはベクトル出力を持つMatlab関数です。ベクトルの各エントリは、変数のいくつかまたは他の評価です。各ステップでMatlabはeqauls 0のうちの1つが0かどうかをチェックします。そのうちの1つが0の場合、DDEが停止し、その部分解を使用してDDEを履歴として再起動する必要があります。

sol = dde23(ddefun, lags, history, [t0 tEnd], options); 

あなたは、ベクトルのエントリの一つは（あなたがxが0に等しいときDDEを停止したいと）xになるだろう。この場合、（予告solとt0変更）

sol = dde23(ddefun, lags, sol, [tCurrent tEnd], options); 

を実行します。また、コードelseif dxTmp <= 0の行は別の不連続を作成するため、dxTmpが0になるイベント、つまり1/(1+y*z) - xがベクター出力の別のコンポーネントになるイベントが必要です。

ODEを再起動すると、Matlabはその時点で不連続性があると自動的に判断します。そのため、Matlabにはそこに存在することを心配する必要はありません。

次の問題は、Matlabのは非常に右のそれを達成したことがない、x、y、zとXの値は若干のマイナスになるということです。問題を作成するために起こっている場合はこのように、あなたは、デリバティブを計算する前に

if x < 0 
    x = 0; 
end 

でxの値（同様に他の値）を補正することになるでしょう。しかし、これはローカルでxの値を変更するだけです。したがって、の最終的なソリューションのxのすべての負の値を0に変更することもできます。 をsol = dde23(ddefun, lags, sol, [tCurrent tEnd], options);に入力する前にそれを変更しようとしないことをお勧めします。私はその試みをいくつか試してみました。

Answer 3

さらにという簡単な答えがあります。createOptimProblemを使用して最適化を設定してください。このメソッドを使用してすべてのパラメータの境界を含める必要がありますが、パラメータが正のままであることを強制するのは簡単です。ここ

詳細：http://www.mathworks.com/help/gads/createoptimproblem.html