Unable to perform assignment because the left and right sides have a different number of elements.

1 visualización (últimos 30 días)
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

Respuesta aceptada

Walter Roberson
Walter Roberson el 9 de Abr. de 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 comentarios
Nana Kwaku Okai-Mensah
Nana Kwaku Okai-Mensah el 9 de Abr. de 2021
I just changed that and that error is gone.
However, now I'm getting this error:
% % Index exceeds the number of array elements (1).
% %
% % Error in LordHlpUs>sir_model (line 55)
% % dS = ((-beta*(B(tt)/(B(tt)+k))*S(tt) + ub*(S(tt)+I(tt)+R(tt)) - ud*S(tt)) * dt)';
[S,I,R,B_] = sir_model(beta,gamma,delta,ub,uc,ud,E,B,k,N,I0,T,dt);
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)
% if delta = 0 we assume a model without immunity loss
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);
B_(1) = B;
for tt = 1:(T/dt)-1
% Equations of the model
dS = ((-beta*(B(tt)/(B(tt)+k))*S(tt) + ub*(S(tt)+I(tt)+R(tt)) - ud*S(tt)) * dt)';
dI = (beta*(B(tt)/(B(tt)+k))*S(tt) - gamma*I(tt) - (uc+ud)*I(tt))* dt;
dR = (gamma*I(tt) - ud*R(tt)) * dt;
dB = (E*I(tt) - delta*B(tt)) * dt;
S(tt+1) = S(tt) + dS;
I(tt+1) = I(tt) + dI;
R(tt+1) = R(tt) + dR;
B_(tt+1) = B_(tt) + dB;
end
end
% what's wrong with my indexing?
Walter Roberson
Walter Roberson el 9 de Abr. de 2021
Your input B is a scalar.
You might be confusing B and B_ (which is a vector)

Iniciar sesión para comentar.

Más respuestas (1)

the cyclist
the cyclist el 9 de Abr. de 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.

Categorías

Más información sobre Statics and Dynamics en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by