Debugging an iterative program
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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.![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/996960/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/996960/image.png)
1 Kommentar
Antworten (1)
Rishi
am 10 Okt. 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1507149/image.jpeg)
If the value of ‘sv(1)’ is checked from the MATLAB Command Window, the output is still 100.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1507154/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1507159/image.jpeg)
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:
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:
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.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Symbolic Math Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!