Getting unexpected results while using ode45

Hi,

I have a system of 10 coupled differential equations. It is a vector-host epidemic model, which captures the transmission of a disease between human population and insect population. Since it is a simple system of differential equations, I am using solvers (ode45) for non-stiff problem type.

There are 10 differential equations, each representing 10 different state variables. There are two functions which have the same system of 10 coupled D.E.’s. One is called “NoEffects_derivative_6_15_2012.m” which contains the original system of D.E.’s. The other function is called “OnlyLethal_derivative_6_15_2012.m” which contains the same system of D.E.’s with an increased withdrawal rate starting at time, gamma=32 days and that withdrawal rate decays exponentially with time.

I use ode45 to solve both the systems, using the same initial conditions. Time vector is also the same for both systems, going from t0 to tfinal. The vector “tspan” contains the time values going from t0 to tfinal, each with a increment of 0.25 days, making a total of 157 time values.

The solution values are stored in matrices “ye0” and ”yeL”. Both these matrices contain 157 rows and 10 columns (for the 10 state variable values).

When I compare the value of the 10th state variable, for the time=tfinal, in the matrix “ye0” and “yeL” by plotting the difference, I find it to be becoming negative for some time values. (using the command: plot(te0,ye0(:,10)-yeL(:,10));). This is not expected. For all time values from t0 till tfinal, the value of the 10 state variable, should be greater, as it is the solution obtained from a system of D.E.’s which did not have an increased withdrawal rate applied to it.

I am told that there is a bug in my matlab code.

I am not sure how to find out that bug. Or maybe the solver in matlab I am using (ode45) is not efficient and does give this kind of problem.

Can anyone help.

I have tried ode23 and ode113 as well, and yet get the same problem.

The figure (2), shows a curve which becomes negative for time values 32 and 34 and this is showing a result which is not expected. This curve should have a positive value throughout, for all time values.

Is there any other forum anyone can suggest ?

Here is the MAIN SCRIPT FILE:

clear memory;
clear all;

global Nc capitalambda muh lambdah del1 del2 p eta alpha1 alpha2 muv lambdav
global dims Q t0 tfinal gamma Ct0 b1 b2 Ct0r b3 H C m_tilda betaHV bitesPERlanding IC
global tspan Hs Cs betaVH k landingARRAY muARRAY 

Nhh=33898857; 
Nvv=2*Nhh; 
Nc=21571585; 
g=354; % number of public health centers in Bihar state
%Fix human parameters
capitalambda= 1547.02;
muh=0.000046142;
lambdah= 0.07; 
del1=0.001331871263014;
del2=0.000288658;  
p=0.24; 
eta=0.0083; 
alpha1=0.044; 
alpha2=0.0217; 
%Fix vector parameters
muv=0.071428; % UNIT:2.13 SANDFLIES DEAD/SAND FLY/MONTH, SOURCE: MUBAYI ET AL., 2010
lambdav=0.05; % UNIT:1.5 TRANSMISSIONS/MONTH, SOURCE: MUBAYI ET AL., 2010

Ct0=0.054;b1=0.0260;b2=0.0610;
Ct0r=0.63;b3=0.0130;

dimsH=6; % AS THERE ARE FIVE HUMAN COMPARTMENTS
dimsV=3; % AS THERE ARE TWO VECTOR COMPARTMENTS
dims=dimsH+dimsV; % THE TOTAL NUMBER OF COMPARTMENTS OR DIFFERENTIAL EQUATIONS

gamma=32; % spraying is done of 1st feb of the year

Q=0.2554;
H=7933615;
C=5392890;

m_tilda=100000; % assumed value 6.5, later I will have to get it for sand flies or mosquitoes
betaHV=66.67/1000000; % estimated value from the short technical report sent by Anuj
bitesPERlanding=lambdah/(m_tilda*betaHV);
betaVH=lambdav/bitesPERlanding;
IC=zeros(dims+1,1); % CREATES A MATRIX WITH DIMS+1 ROWS AND 1 COLUMN WITH ALL ELEMENTS AS ZEROES

t0=1;
tfinal=40;
for j=t0:1:(tfinal*4-4)
    tspan(1)= t0;
    tspan(j+1)= tspan(j)+0.25;
end
clear j;

% INITIAL CONDITION OF HUMAN COMPARTMENTS
q1=0.8; q2=0.02; q3=0.0005; q4=0.0015;
IC(1,1) = q1*Nhh;
IC(2,1) = q2*Nhh;
IC(3,1) = q3*Nhh;
IC(4,1) = q4*Nhh;
IC(5,1) = (1-q1-q2-q3-q4)*Nhh;
IC(6,1) = Nhh;
% INTIAL CONDITIONS OF THE VECTOR COMPARTMENTS
IC(7,1) = 0.95*Nvv; %80 PERCENT OF TOTAL ARE ASSUMED AS SUSCEPTIBLE VECTORS
IC(8,1) = 0.05*Nvv; %20 PRECENT OF TOTAL ARE ASSUMED AS INFECTED VECTORS
IC(9,1) = Nvv;
IC(10,1)=0;

Hs=2000000;
Cs=3000000;
k=1;
landingARRAY=zeros(tfinal*50,2);
muARRAY=zeros(tfinal*50,2);

[te0 ye0]=ode45(@NoEffects_derivative_6_15_2012,tspan,IC);
[teL yeL]=ode45(@OnlyLethal_derivative_6_15_2012,tspan,IC);

figure(1)
subplot(4,3,1); plot(te0,ye0(:,1),'b-',teL,yeL(:,1),'r-');
xlabel('time'); ylabel('S');
legend('susceptible humans');
subplot(4,3,2); plot(te0,ye0(:,2),'b-',teL,yeL(:,2),'r-');
xlabel('time'); ylabel('I');
legend('Infectious Cases');
subplot(4,3,3); plot(te0,ye0(:,3),'b-',teL,yeL(:,3),'r-');
xlabel('time'); ylabel('G');
legend('Cases in Govt. Clinics');
subplot(4,3,4); plot(te0,ye0(:,4),'b-',teL,yeL(:,4),'r-');
xlabel('time'); ylabel('T');
legend('Cases in Private Clinics');
subplot(4,3,5); plot(te0,ye0(:,5),'b-',teL,yeL(:,5),'r-');
xlabel('time'); ylabel('R');
legend('Recovered Cases');

subplot(4,3,6);plot(te0,ye0(:,6),'b-',teL,yeL(:,6),'r-');
hold on; 
plot(teL,capitalambda/muh); 
xlabel('time'); ylabel('Nh');
legend('Nh versus time');hold off;

subplot(4,3,7); plot(te0,ye0(:,7),'b-',teL,yeL(:,7),'r-');
xlabel('time'); ylabel('X');
legend('Susceptible Vectors');

subplot(4,3,8); plot(te0,ye0(:,8),'b-',teL,yeL(:,8),'r-');
xlabel('time'); ylabel('Z');
legend('Infected Vectors');

subplot(4,3,9); plot(te0,ye0(:,9),'b-',teL,yeL(:,9),'r-');
xlabel('time'); ylabel('Nv');
legend('Nv versus time');

subplot(4,3,10);plot(te0,ye0(:,10),'b-',teL,yeL(:,10),'r-');
xlabel('time'); ylabel('FS');
legend('Total number of human infections');

figure(2)
plot(te0,ye0(:,10)-yeL(:,10));
xlabel('time'); ylabel('FS(without intervention)-FS(with lethal effect)');
legend('Diff. bet. VL cases with and w/o intervention:ode45');

THE FUNCTION FILE: NoEffects_derivative_6_15_2012
function dx=NoEffects_derivative_6_15_2012(t,x)

global Nc capitalambda muh del1 del2 p eta alpha1 alpha2 muv    
global dims m_tilda betaHV bitesPERlanding betaVH     

dx=zeros(dims+1,1);
% t
% dx
dx(1,1)= capitalambda-(m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-muh*x(1,1);
dx(2,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-(del1+eta+muh)*x(2,1);
dx(3,1)=p*eta*x(2,1)-(del2+alpha1+muh)*x(3,1);
dx(4,1)=(1-p)*eta*x(2,1)-(del2+alpha2+muh)*x(4,1);
dx(5,1)=alpha1*x(3,1)+alpha2*x(4,1)-muh*x(5,1);
dx(6,1)=capitalambda -del1*x(2,1)-del2*x(3,1)-del2*x(4,1)-muh*x(6,1);
dx(7,1)=muv*(x(7,1)+x(8,1))-bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-muv*x(7,1);
% dx(8,1)= lambdav*x(7,1)*x(2,1)/(x(6,1)+Nc)-muvIOFt(t)*x(8,1);
dx(8,1)= bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-muv*x(8,1);
dx(9,1)= (muv-muv)*x(9,1);
dx(10,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/x(9,1);

THE FUNCTION FILE: OnlyLethal_derivative_6_15_2012
function dx=OnlyLethal_derivative_6_15_2012(t,x)
global Nc capitalambda muh del1 del2 p eta alpha1 alpha2 muv    
global dims m_tilda  betaHV bitesPERlanding betaVH k muARRAY

dx=zeros(dims+1,1);

% the below code saves some values into the second column of the two arrays
% t
muARRAY(k,1)=t;
muARRAY(k,2)=artificialdeathrate1(t);
k=k+1;
dx(1,1)= capitalambda-(m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-muh*x(1,1);

dx(2,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/(x(7,1)+x(8,1))-(del1+eta+muh)*x(2,1);

dx(3,1)=p*eta*x(2,1)-(del2+alpha1+muh)*x(3,1);

dx(4,1)=(1-p)*eta*x(2,1)-(del2+alpha2+muh)*x(4,1);

dx(5,1)=alpha1*x(3,1)+alpha2*x(4,1)-muh*x(5,1);

dx(6,1)=capitalambda -del1*x(2,1)-del2*( x(3,1)+x(4,1) ) - muh*x(6,1);

dx(7,1)=muv*( x(7,1)+x(8,1) )- bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc) - (artificialdeathrate1(t) + muv)*x(7,1);

dx(8,1)= bitesPERlanding*betaVH*x(7,1)*x(2,1)/(x(6,1)+Nc)-(artificialdeathrate1(t) + muv)*x(8,1);

dx(9,1)= -artificialdeathrate1(t) * x(9,1);

dx(10,1)= (m_tilda)*bitesPERlanding*betaHV*x(1,1)*x(8,1)/x(9,1);

THE FUNCTION FILE: artificialdeathrate1
function art1=artificialdeathrate1(t)
global Q Hs H Cs C    

art1= Q*Hs*iOFt(t)/H + (1-Q)*Cs*oOFt(t)/C ;

THE FUNCTION FILE: iOFt(t)
function i=iOFt(t)
global gamma tfinal Ct0 b1 

if (t>=gamma && t<=tfinal)
  i= Ct0*exp(-b1*(t-gamma));
else
  i=0;
end

THE FUNCTION FILE: oOFt(t)
function o=oOFt(t)

global gamma Ct0 b2 tfinal

if (t>=gamma && t<=(tfinal))
  o = Ct0*exp(-b2*(t-gamma));
else
  o = 0;
end