Unable to perform assignment because the left and right sides have a different number of elements.
2 views (last 30 days)
Show older comments
Nana Kwaku Okai-Mensah
on 9 Apr 2021
Commented: Walter Roberson
on 9 Apr 2021
Please help me out.
% Model parameters
beta = 0.07; % daily rate of infection
gamma = 0.12; % rate of recovery
delta = 0.02; % rate of immunity loss
N = 15*10^3; % Total population N = S + I + R
I0 = 10; % initial number of infected
T = 100; % period of 300 days
B = 1000; %
k = 10^5; %
ub = 63*10^-6; % daily rate of birthsper thousand
uc = 68.49*10^-6; % rate of deaths related to cholera
ud = 24.66*10^-6; % rate of deaths unrelated to cholrea
E = 10.0; % cells per mL water per day
dt = 1/4; % time interval of 6 hours (1/4 of a day)
fprintf('Value of parameter R0 is %.2f',N*beta/gamma)
% Calculate the model
[S,I,R,B_] = sirb_model(beta,gamma,delta,ub,uc,ud,E,B,k,N,I0,T,dt);
% Plots that display the epidemic outbreak
tt = 0:dt:T-dt;
% Curve
plot(tt,S,'b',tt,I,'r',tt,R,'g','LineWidth',2); grid on;
xlabel('Days'); ylabel('Number of individuals');
legend('S','I','R');
plot(tt,B,'b','LineWidth',2); grid on;
xlabel('Days'); ylabel('Number of Viruses');
legend('B');
% Map
plot(I(1:(T/dt)-1),I(2:T/dt),"LineWidth",1,"Color",'r');
hold on; grid on;
plot(I(2),I(1),'ob','MarkerSize',4);
xlabel('Infected at time t'); ylabel('Infected at time t+1');
hold off;
Function that calculates the model evolution over a time period T.
function [S,I,R,B_] = sir_model(beta,gamma,delta,ub,uc,ud,E,B,k,N,I0,T,dt)
S = zeros(1,T/dt);
S(1) = N;
I = zeros(1,T/dt);
I(1) = I0;
R = zeros(1,T/dt);
B_ = zeros(1,T/dt);
for tt = 1:(T/dt)-1
% Equations of the model
dS = ((-beta*(B/(B+k))*S(tt) + ub*(S(tt)+I(tt)+R(tt)) - ud*S) * dt)';
dI = (beta*(B/(B+k))*S(tt) - gamma*I(tt) - (uc+ud))* dt;
dR = (gamma*I(tt) - ud*R(tt)) * dt;
dB = (E*I - delta*B) * dt;
S(tt+1) = S(tt) + dS; %%%% The error occurs around here
I(tt+1) = I(tt) + dI;
R(tt+1) = R(tt) + dR;
B_(tt+1) = B_(tt) + dB;
end
end
0 Comments
Accepted Answer
Walter Roberson
on 9 Apr 2021
S = zeros(1,T/dt);
S(1) = N;
S is a vector.
dS = ((-beta*(B/(B+k))*S(tt) + ub*(S(tt)+I(tt)+R(tt)) - ud*S) * dt)';
^
You are using all of S in the marked place, so dS is not a scalar.
I = zeros(1,T/dt);
I(1) = I0;
I is a vector.
dB = (E*I - delta*B) * dt;
^
You are using all of I in the marked place, so dB is not a scalar.
2 Comments
Walter Roberson
on 9 Apr 2021
Your input B is a scalar.
You might be confusing B and B_ (which is a vector)
More Answers (1)
the cyclist
on 9 Apr 2021
It looks like inside your function, in these lines
dS = ((-beta*(B/(B+k))*S(tt) + ub*(S(tt)+I(tt)+R(tt)) - ud*S) * dt)';
dI = (beta*(B/(B+k))*S(tt) - gamma*I(tt) - (uc+ud))* dt;
dR = (gamma*I(tt) - ud*R(tt)) * dt;
dB = (E*I - delta*B) * dt;
you probably intended to use S(tt) in the first line, and I(tt) in the last line, instead of just S and I.
See Also
Categories
Find more on Biological and Health Sciences in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!