非負の ODE の解

このトピックでは、ODE の解が非負になるように制約する方法を説明します。非負の制約を課すことは必ずしも簡単ではありませんが、方程式の物理的な解釈や解の性質によっては必要になる場合があります。解に対するこのような制約は、必要な場合にのみ行ってください。たとえば、この制約がないと積分が失敗する場合や解が適切ではない場合です。

解の特定の要素が非負でなければならない場合、odeset を使用して、これらの要素のインデックスに NonNegative オプションを設定します。このオプションは、ode23s または ode15i には使用できません。また、質量行列を使用する問題に対して適用する陰的ソルバー (ode15s、ode23t、ode23tb) にも使用できません。特に、必ず特異質量行列をもつ DAE 問題に非負の制約を課すことはできません。

例: 絶対値関数

次の初期値問題について考えます。

$$y^{\prime }=-|y|,$$

この問題では、初期条件を $$y(0)=1$$ として、区間 $$[0,40]$$ について解を求めます。この ODE の解はゼロに減衰します。ソルバーが解として負の値を生成する場合は、この値を使用して ODE の解の追跡が開始され、計算解が $$-\infty$$ に発散すると、計算が最終的に失敗します。NonNegative オプションを使用すると、このような積分の失敗を防ぐことができます。

ODE の求解で ode45 を追加オプションなしで使用した場合と、これを NonNegative オプションを設定した場合について、$$y(t)=e^{-t}$$ の解析解との比較を行います。

ode = @(t,y) -abs(y);

% Standard solution with |ode45|
options1 = odeset('Refine',1);
[t0,y0] = ode45(ode,[0 40],1,options1);

% Solution with nonnegative constraint
options2 = odeset(options1,'NonNegative',1);
[t1,y1] = ode45(ode,[0 40],1,options2);

% Analytic solution
t = linspace(0,40,1000);
y = exp(-t);

この 3 つの解をプロットして比較します。解が $$-\infty$$ に向かわないようにするには、非負の制約を課すことが非常に重要です。

plot(t,y,'b-',t0,y0,'ro',t1,y1,'k*');
legend('Exact solution','No constraints','Nonnegativity', ...
       'Location','SouthWest')

例: Knee 問題

非負の解が必要な問題のもう一つの例として、"Knee 問題" を取り上げます。この問題はサンプル ファイル kneeode にコード化されています。この方程式は次のとおりです。

$$ϵy^{\prime }=(1-x)y-y^2,$$

この問題では、初期条件を $$y(0)=1$$ として、区間 $$[0,2]$$ について解を求めます。パラメーター $$ϵ$$ は通常、$$0<ϵ\ll 1$$ を満たすものとされ、この問題では $$ϵ=1\times 10^{-6}$$ が使用されます。この ODE の解は、$$x<1$$ の場合は $$y=1-x$$ に近付き、$$x>1$$ の場合は $$y=0$$ に近付きます。ただし、既定の許容誤差を使用して数値解を計算すると、この解は積分の全区間に対して $$y=1-x$$ のアイソクラインをなすことが示されます。非負という制約を課すと正確な解になります。

非負の制約を課す場合と課さない場合の Knee 問題の解を求めます。

epsilon = 1e-6;
y0 = 1;
xspan = [0 2];
odefcn = @(x,y,epsilon) ((1-x)*y - y^2)/epsilon;

% Solve without imposing constraints
[x1,y1] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0);

% Impose a nonnegativity constraint
options = odeset('NonNegative',1);
[x2,y2] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0, options);

これらの解をプロットして比較します。

plot(x1,y1,'ro',x2,y2,'b*')
axis([0,2,-1,1])
title('The "knee problem"')
legend('No constraints','Non-negativity')
xlabel('x')
ylabel('y')

参照

[1] Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne, "Non-negative solutions of ODEs," Applied Mathematics and Computation Vol. 170, 2005, pp. 556-569.

参考

odeset

関連するトピック

• ODE ソルバーの選択
• ODE オプションの概要