Debugging an iterative program

Hi,

Iam using an iterative method to solve a SEIRV mathematical model.

all is well but, when the compiler start the for loop which, I have used for loop to iterate time and other matrices.

The first loop alone runs successfully but, from the second loop the previous values are zeroes which is not the expected results

Any help is appreciated.

clc

close all

clear all

%%

%%Initial values

t0 = 0;

tend = 50;

%tspan=[t0,tend];

y0 = [100,10,3,0,0];

%%

%%Constant(rates) values

A=8859.23*(10^(4)); % Influx of population

be=6.11*(10^(-8)); % Contact rate between E and S

bi=2.62*(10^(-8)); % Contact rate between I and S

bv=3.03*(10^(-8)); % Transmission from envi to humans

m=3.01*(10^(-2)); % Natural death population

a=0.143; % quarantine period of infected

o=0.01; % Death rate from Infected

g=0.67; % Recovery rate recovered from disease and alive

p1=1.3; % Exposed population which contributing the virus in the envi

p2=0.06; % Infected population which contributing the virus in the envi

th=2.0; % Diesease-induced removal rate

%%

h= 1;

tf= (tend/h);

al = 0.83; % Fractional order parameter

w = 1; % Normalization Function

s = zeros(1,tf);

e = zeros(1,tf);

i = zeros(1,tf);

r = zeros(1,tf);

v = zeros(1,tf);

t = zeros(tf,1); % time (row vector)

sv = zeros(tf,5);

for ii = 1:tf

t(ii,1) = ( t0 + ((ii-1)*h));

if ii == 1 %% Assigning the initial values to the sv matrix where, sv(:,1)= S ; sv(:,2) = E, sv (:,3) = I; sv(:,4) = R, sv(:,5) = V

sv(ii,1) = y0(1);

sv(ii,2) = y0(2);

sv(ii,3) = y0(3);

sv(ii,4) = y0(4);

sv(ii,5) = y0(5);

s(1,ii) = A-(be*sv(ii,1)*sv(ii,2))-(bi*sv(ii,1)*sv(ii,3))- (bv*sv(ii,1)*sv(ii,5))- (m*sv(ii,1)); %% Assigning the first calculated rates of S,E,I,R,V into the respected rate vector.

e(1, ii) = (be*sv(ii,1)*sv(ii,2))+(bi*sv(ii,1)*sv(ii,3))+ (bv*sv(ii,1)*sv(ii,5)) - ((a+m)*sv(ii,2));

i(1,ii) = (a*sv(ii,2))-((o+m+g)*sv(ii,3));

r(1,ii) = (g*sv(ii,3))-(m*sv(ii,4));

v (1,ii) = (p1*sv(ii,2))+p2-(th*sv(ii,5));

elseif ii == 2

sv(ii,1) = sv((ii-1),1) + (((1-al)/w)+((3*al*h)/(2*w)))*s(1,(ii-1)); %% Assigning the second element in the respective columns (S,E,I,R,V)

sv(ii,2) = sv((ii-1),2) + (((1-al)/w)+((3*al*h)/(2*w)))*e(1,(ii-1)); %% As only values @tn is available, the term involving the values @t(n-1) is neglected

sv(ii,3) = sv((ii-1),3) + (((1-al)/w)+((3*al*h)/(2*w)))*i(1,(ii-1));

sv(ii,4) = sv((ii-1),4) + (((1-al)/w)+((3*al*h)/(2*w)))*r(1,(ii-1));

sv(ii,5) = sv((ii-1),5) + (((1-al)/w)+((3*al*h)/(2*w)))*v(1,(ii-1));

s(1,ii) = A-(be*sv(ii,1)*sv(ii,2))-(bi*sv(ii,1)*sv(ii,3))- (bv*sv(ii,1)*sv(ii,5))- (m*sv(ii,1)); %% Assigning the calculated rates of S,E,I,R,V into the respected rate vector

e(1, ii) = (be*sv(ii,1)*sv(ii,2))+(bi*sv(ii,1)*sv(ii,3))+ (bv*sv(ii,1)*sv(ii,5)) - ((a+m)*sv(ii,2));

i(1,ii) = (a*sv(ii,2))-((o+m+g)*sv(ii,3));

r(1,ii) = (g*sv(ii,3))-(m*sv(ii,4));

v (1,ii) = (p1*sv(ii,2))+p2-(th*sv(ii,5));

else

sv(ii,1) = sv((ii-1),1) + (((1-al)/w)+((3*al*h)/(2*w)))*s(1,(ii-1)) + (((1-al)/w)+((al*h)/(2*w)))*s(1,(ii-2)); %% Assigning the every other element in the respective columns (S,E,I,R,V)

sv(ii,2) = sv((ii-1),2) + (((1-al)/w)+((3*al*h)/(2*w)))*e(1,(ii-1)) + (((1-al)/w)+((al*h)/(2*w)))*e(1,(ii-2));

sv(ii,3) = sv((ii-1),3) + (((1-al)/w)+((3*al*h)/(2*w)))*i(1,(ii-1)) + (((1-al)/w)+((al*h)/(2*w)))*i(1,(ii-2)); %% The whole equation is employed here as values @tn and @t(n-1) is known

sv(ii,4) = sv((ii-1),4) + (((1-al)/w)+((3*al*h)/(2*w)))*r(1,(ii-1)) + (((1-al)/w)+((al*h)/(2*w)))*r(1,(ii-2));

sv(ii,5) = sv((ii-1),5) + (((1-al)/w)+((3*al*h)/(2*w)))*v(1,(ii-1)) + (((1-al)/w)+((al*h)/(2*w)))*v(1,(ii-2));

s(1,ii) = A-(be*sv(ii,1)*sv(ii,2))-(bi*sv(ii,1)*sv(ii,3))- (bv*sv(ii,1)*sv(ii,5))- (m*sv(ii,1));

e(1, ii) = (be*sv(ii,1)*sv(ii,2))+(bi*sv(ii,1)*sv(ii,3))+ (bv*sv(ii,1)*sv(ii,5)) - ((a+m)*sv(ii,2)); %% Assigning the calculated rates of S,E,I,R,V into the respected rate vector

i(1,ii) = (a*sv(ii,2))-((o+m+g)*sv(ii,3));

r(1,ii) = (g*sv(ii,3))-(m*sv(ii,4));

v (1,ii) = (p1*sv(ii,2))+p2-(th*sv(ii,5));

end

end

%% Plotting time vs. S,E,I,R,V

figure(1)

clf

plot(t(:,1),sv(:,1),'LineWidth', 3); title('Susceptible Population vs Time')

xlabel('Time')

ylabel('Susceptible population')

figure(2)

clf

plot(t,sv(:,2),'LineWidth', 3); title('Exposed Population vs Time')

xlabel('Time')

ylabel('Exposed population')

figure(3)

clf

plot(t,sv(:,3),'LineWidth', 3); title('Infected Population vs Time')

xlabel('Time')

ylabel('Infected population')

figure(4)

clf

plot(t,sv(:,4),'LineWidth', 3); title('Recovered Population vs Time')

xlabel('Time')

ylabel('Recovered population')

figure(5)

clf

plot(t,sv(:,5),'LineWidth', 3); title('Concentration vs Time')

xlabel('Time')

ylabel('Death population')

This is the results for the first loop,first two loop and 1 to tf loops. in the first pic the values are assigned in the first rows but from the second loop the first row becomes zeroes but the value from the first loop should be there in the first row. Then in the last pic after a certain point the values are registered as NaN.

[Bjorn Gustavsson]

What is the biological justification for a constant recovery-rate?

Answer 1

[Rishi]

Hi Kesavan,

I understand that you want to know why the first row of the matrix ‘sv’ is changing in the MATLAB Command Window after the second iteration, and why the last few rows show up as 'NaN'.

After running the code, when you check the variable ‘sv’ through the MATLAB Workspace, you can see that the values are unchanged. Also, the values have the order of 10^186 from the 14th row. When displayed in the MATLAB Command Window, the entire matrix ‘sv’ is multiplied by 1.0e+186, and the values in the initial rows become too small to be displayed in three digits after the decimal. If the value of ‘sv(1)’ is checked from the MATLAB Command Window, the output is still 100.

If the value of ‘sv(1)’ is checked from the MATLAB Command Window, the output is still 100.

The range for a negative number of type double is between -1.79769 x 10308 and -2.22507 x 10-308, and the range for positive numbers is between 2.22507 x 10-308 and 1.79769 x 10308, as mentioned in the following MATLAB documentation:

https://www.mathworks.com/help/releases/R2021b/matlab/ref/double.html?s_tid=doc_ta#:~:text=The%20range%20for%20a%20negative%20number%20of%20type%20double%20is%20between%20%2D1.79769%20x%2010308%20and%20%2D2.22507%20x%2010%2D308%2C%20and%20the%20range%20for%20positive%20numbers%20is%20between%202.22507%20x%2010%2D308%20and%201.79769%20x%2010308.

Any further operations on the matrix ‘sv’ result in values which are out of this range, and hence, the answer comes out to be ‘inf’. Then, while calculating the values of 's', results in thesubtraction of ‘inf’ from ‘inf’, giving the output as ‘NaN’.

To work with values beyond the range of double, symbolic variables can be used. You can learn more about them in the following documentation:

https://www.mathworks.com/help/releases/R2021b/symbolic/sym.html

Here is the edited initial section of code that generates the desired result:

clc
close all
clear all

%% Initial values
t0 = 0;
tend = 50;
y0 = sym([100, 10, 3, 0, 0]);
h = 1;
tf1 = tend/h;
%% Constant (rate) values
A = sym(8859.23 * (10^4));  % Influx of population
be = sym(6.11 * (10^(-8)));   % Contact rate between E and S
bi = sym(2.62 * (10^(-8)));   % Contact rate between I and S
bv = sym(3.03 * (10^(-8)));   % Transmission from environment to humans
m = sym(3.01 * (10^(-2)));    % Natural death population
a = sym(0.143);             % Quarantine period of infected
o = sym(0.01);              % Death rate from Infected
g = sym(0.67);              % Recovery rate recovered from disease and alive
p1 = sym(1.3);              % Exposed population contributing the virus in the environment
p2 = sym(0.06);             % Infected population contributing the virus in the environment
th = sym(2.0);                % Diesease-induced removal rate
%%    
h = sym(1);
al = sym(0.83);         % Fractional order parameter
w = sym(1);             % Normalization Function


s = sym(zeros(1, tf1));
e = sym(zeros(1, tf1));
i = sym(zeros(1, tf1));
r = sym(zeros(1, tf1));
v = sym(zeros(1, tf1));

t = sym(zeros(tf1, 1));   % time (row vector)
sv = sym(zeros(tf1, 5));

Hope this helps.