Making time steps discrete for ode45

Hi all,

I have a code in which I am solving a single Differential equation by using ode45. I have time going from t0=1 day to tfinal=365 days. The ode45 solves a DE that has vector, which has 365 values in it. So, I make tspan have increments of one day before I pass it to ode45. Even then I noticed that ode45 tries to solve the DE for continuous time values. I placed the variable name, "t" before the DE in the function, named single_derivative.m to see for which values of t does ode45 solve the DE. I want it to solve the DE for only the 365 days. The main script file and the function file are given below:

clear memory
clear all

global Q tfinal epsilon gamma Ct0 b1 b2 Ct0r b3 muv percent_landing
global H C Hs Cs Nv h_bar c_bar  

% generating the values of exponential decay of lethal effect indoors
tfinal=365; % user gives the number of days for which the simulation will be run
epsilon=90; % user gives the number of days for which the insecticide effect remains

gamma=32; % spraying is done of 1st feb of the year
Ct0=0.054;b1=-0.0260;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        h(t)= Ct0*exp(b1*tau);
    else
        h(t)=0;
    end
end

% generating the values of exponential decay of lethal effect outdoors
Ct0=0.054;b2=-0.0610;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        z(t)= Ct0*exp(b2*tau);
    else
        z(t)=0;
    end
end

% generating the values of exponential decay of repellent effect 
Ct0r=0.63;b3=-0.0130;
for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
       tau=t-gamma;
        r(t)= Ct0*exp(b3*tau);
    else
        r(t)=0;
    end
end

% net vector death rate with intervention
muv=0.071;
percent_landing=0.65;
Q=0.2554;
H=7933615;
C=5392890;
Hs=3000000;
Cs=2000000;

for t=1:tfinal
    if(t>gamma && t<(gamma+epsilon))
        tau=t-gamma;
        muvI(t)= Q*(Hs/H)*(1-percent_landing)*r(t)*muv + Q*(Hs/H)*(1-percent_landing)*(1-r(t))*(h(t)+muv)+ Q*(Hs/H)*(percent_landing)*(h(t)+muv)+ Q*(1-Hs/H)*muv + (1-Q)*(Cs/C)*(1-percent_landing)*r(t)*muv + (1-Q)*(Cs/C)*(1-percent_landing)*(1-r(t))*(z(t)+muv)+(1-Q)* (Cs/C)*(percent_landing)* (z(t)+muv)+(1-Q)*(1-Cs/C)*muv;
    else 
        muvI(t)= muv;
    end
end 

% Adjusted sand fly population gamma days after insecticide application
% Note: the values are generated for discrete days of the years, denoted by n 
Nv=10000; % assumed the vector population as constant before spray
for n=1:tfinal
    if(n>=1 && n<=(gamma-1))
        Nv_adjusted(n)=Nv;
    else
        Nv_adjusted(n)=Nv_adjusted(n-1)*(1-muvI(n)+muv);
    end 
end

plot(h)
plot(z)
plot(r)
plot(muvI)
plot(Nv_adjusted)

% rate of change of adjusted vector population
for t=1:tfinal
    if (t>=gamma && t<=tfinal)
        del_Nv_adjusted(t)=(muvI(t)-muv)*Nv_adjusted(t);
    else
        del_Nv_adjusted(t)=0;
    end
end

% rate of change of landings averted
h_bar=5.4;
c_bar=4;

for t=1:tfinal
    if (t>=gamma && t<=tfinal)
        del_landings_averted(t)=Q*del_Nv_adjusted(t)*(Hs/H)*(1-percent_landing)*(r(t))*(t-gamma)/(Hs*h_bar)+ (1-Q)*del_Nv_adjusted(t)*(Cs/C)*(1-percent_landing)*(r(t))*(t-gamma)/(Cs*c_bar);
    else
        del_landings_averted(t)=0;
    end
end

plot(del_Nv_adjusted)
plot(del_landings_averted)

% Need to write code for landings_averted, L(t) using integration
t0=1;
tfinal=374;
tspan=[t0:1:tfinal];
ICNv_adjusted=zeros(1,1); 
ICNv_adjusted(1,1) = 5000; % arbitrarily taken later will take the value of the solution of ODE45 for time equal to gamma

[t y]=ode45(@single_derivative, tspan, ICNv_adjusted, [], muvI, muv);
% end of main script file

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Function file named: single_derivative.m

function dx=single_derivative(t,x,muvI,muh)

t

dx=zeros(1,1);

dx(t,1) = (muh-muvI(t))*x(1,1); 
% end of function file

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

I am getting the following error on the command prompt:

t =

       1

t =

     1.200000000000000

??? Attempted to access muvI(1.2); index must be a positive integer or logical.

Error in ==> single_derivative at 7
dx(t,1) = (muh-muvI(t))*x(1,1);
Error in ==> ode45 at 324
    f(:,2) = feval(odeFcn,t+hA(1),y+f*hB(:,1),odeArgs{:});

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

ode45 is trying to solve the DE for t=1 and then t=1.2 and actually I there is not value for the vector muvI(t) for t=1.2

This vector has only 365 values. So I want ode45 to solve and give solutions of the DE for only discrete time values from 1 to 365.

Please suggest on how to solve this problem. Is it appropriate to make changes in the ode45 function, which is a built in function? Or should I use another solver.

Thanks in advance.

Kaushik.

Answer 1

[Jan]

ODE45 solves continuos functions only. Although you could cheat using something like (muh-muvI(round(t))) * x(1,1), the quality of the result will suffer from the discontinuities.

In addition the left hand side of dx(t,1) = (muh-muvI(t))*x(1,1) looks strange. It looks like the dimension of the state variable and the time steps are confused. In the line before you define dx as a scalar.

A standard method to solve discontinuos functions with ODE45 is to restart the integration at each step.