0

私はMatLabでいくつかのODEを解こうとしていますが、方程式の変数は母集団であるため、正になるように制約する必要があります。そこで、方程式ソルバーを呼び出して非負にする前に odeset() を使用してみましたが、後で値をプロットすると、実際には負になることがあります (以下のコードでは、マゼンタの線です)。私は何を間違っていますか？

ここにいくつかのコードがあります：

%Lots of variables
includeJ=1;
cullLIRate=1/2000;
cullDIRate=1/2000;
N = 16800;
beta = 2e-7;
delta = 0.5;
gamma = 1/50;
sigma = 1/400;
mu = 1/365;
maxTime = 30*365;
kappa = N;
gR = 0.05;
mJ = 1/3650;
initJPerAdult = 10;
numInitE = 1000;
TSpan = [0,maxTime];

initState = [N-numInitE,numInitE,0,0,0,initJPerAdult*N];

%IMPORTANT BIT HERE
options = odeset('NonNegative', 1:6)
scirSoln = ode45(@equation,TSpan,initState,[],beta,delta,gamma,sigma,mu,kappa,gR,mJ,cullLIRate,cullDIRate,includeJ);

scirVals = deval(scirSoln,timeToPlot);
plot(timeToPlot,scirVals(1,:));
hold on;
plot(timeToPlot,scirVals(3,:),'k');
plot(timeToPlot,scirVals(4,:),'g');
plot(timeToPlot,scirVals(6,:),'m');

timeToPlot = [0:max(TSpan)/1000:max(TSpan)];

方程式(...) のコードは次のとおりです。

function retVal = equation(t,y,beta,delta,gamma,sigma,mu,kappa,gR,mJ,cullLIRate,cullDIRate,includeJ)
retVal = zeros(6,1);

S = y(1);
E = y(2);
LI = y(3);
DI = y(4);
R = y(5);
J = y(6);

retVal(1)= mJ * J - beta * S * (delta * LI + DI);
retVal(2) = beta * S * (delta * LI + DI) - gamma * E;
retVal(3) = gamma * E - (cullLIRate + sigma) * LI;
retVal(4) = sigma * LI - (mu + cullDIRate) * DI;
retVal(5) = mu * DI + cullLIRate* LI + cullDIRate * DI;
retVal(6) = gR * S * (1 - S / kappa) - mJ * J;

end
4

1 に答える 1

1

定義したodeset（オプション変数）をODE45-ソルバーに渡していない。

ODE45の構文は次のとおりです。[T、Y] = ODE45（ODEFUN、TSPAN、Y0、OPTIONS、P1、P2 ...）

うまくいったことをうれしく思います！:)

于 2013-02-24T18:32:58.780 に答える