James Elderfield
James Elderfield

Reputation: 2507

Using NonNegative setting in Matlab odeset()

I'm trying to solve some ODEs in MatLab and seeing as the variables in the equations are populations they need to be constrained to being positive. So I tried using odeset() before calling the equation solver to make them non-negative but on plotting the values afterwards they are actually negative at times (in the code below it is the magenta line). What am I doing wrong?

Here's some code:

%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)];

The code for equation(...) is:

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

Upvotes: 0

Views: 4741

Answers (1)

Groot
Groot

Reputation: 14251

You are not passing your defined odeset (options variable) to the ODE45 - solver.

The syntax for the ODE45 is: [T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2...)

Glad it worked! :)

Upvotes: 1

Related Questions