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

The code shown below is my attempt to develop an SEIR model to solve the system of Ordinary Differential Equations using RK4. An error pops up when running the code stating :

"Unable to perform assignment because the left and right sides have a

different number of elements.

Error in SEIR_RK4 (line 55)

S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6;"

I've looked through numerous sources but I'm still puzzeled how to make the code run.

r0 = 3;
R = r0;%For Now
Tinc = 0.5; %Time of Incubation ~ 2 days
Tinf = 0.1; %Time of Infected ~ 10 days
t0 = 0; %Initial Time - Given
tf = 300;%Final Time - Given
N = 300; %Number of Points (or "Day Intervals") % M = N 
tol = 10^(-6); %Tolerance

h = (tf-t0)/N;%Step Size
t = [t0:h:tf]';%Discretized time step

S0 = 0.8; %Initial amount of Susceptibles
E0 = 0.2; %Initial amount of Incubants
I0 = 0.0; %Initial amount of Infectants
R0 = 0.0; %Initial amount of Recovered

%Initial Values X(0) - Matrix (4x1)
S = zeros(1,N); S(1) = S0;
E = zeros(1,N); E(1) = E0;
I = zeros(1,N); I(1) = I0;
R = zeros(1,N); R(1) = R0;

% f(t,x) = dx/dt
%f = @(t,x) [(-R*x(3)*x(1)/Tinf);(R*x(3)*x(1)/Tinf)-(x(2)/Tinc);(x(2)/Tinc)-(x(3)/Tinf);(x(3)/Tinf)];
dS = @(s,e,i,r) (-R.*i.*s/Tinf);
dE = @(s,e,i,r) (R.*i.*s/Tinf)-(e/Tinc);
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);
dR = @(s,e,i,r) (i/Tinf);

for i = 1:N
    %k1 = h*f(tn,yn)
    s1 = h*dS( S(i), E(i), I(i), R(i));
    e1 = h*dE( S(i), E(i), I(i), R(i));
    i1 = h*dI( S(i), E(i), I(i), R(i));
    r1 = h*dR( S(i), E(i), I(i), R(i));
    %k2 = h*f(tn + (h/2),yn + k1.*(h/2))
    s2 = h*dS( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    e2 = h*dE( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    i2 = h*dI( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    r2 = h*dR( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    %k3 = h*f(tn + (h/2),yn + k2.*(h/2))
    s3 = h*dS( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    e3 = h*dE( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    i3 = h*dI( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    r3 = h*dR( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    %k4 = h*f(tn + h,yn + h.*k3)
    s4 = h*dS( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    e4 = h*dE( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    i4 = h*dI( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    r4 = h*dR( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    %yn+1 = yn + (1/6)*h*(k1+2*k2+2*k3+k4)
    S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6;
    E(i+1) = E(i) + (e1 + 2*e2 + 2*e3 + e4)/6;
    I(i+1) = I(i) + (i1 + 2*i2 + 2*i3 + i4)/6;
    R(i+1) = R(i) + (r1 + 2*r2 + 2*r3 + r4)/6;
end

Answer 1

[Walter Roberson]

Okay, so it turns out that you do not need to use cells. What you do need is to fix

dS = @(s,e,i,r) (-R.*i.*s/Tinf);

R is not

R = r0;%For Now

like you claim, because you later have

R = zeros(1,N); R(1) = R0;

so R is the vector for keeping track of the number of recovered individuals.

You need to rename one of the uses of R, as either the transmittability or the number of recovered individuals.

If you define

odefun = @(t,x) [dS(x(1),x(2),x(3),x(4)); dE(x(1),x(2),x(3),x(4)); dI(x(1),x(2),x(3),x(4)); dR(x(1),x(2),x(3),x(4))];
x0 = [S0; E0; I0; R0];
tspan = t;
[t, x] = ode23(odefun, tspan, x0);

then the results you get would suggest that your implementation of ode23 is incorrect. But also as you look at the E column (second column) and I column (third column) since the values cycle between negative and positive, that suggests that your dE and dI are incorrect.

Question: are your values intended to be proportions or intended to be absolute numbers? You do not have a population size, so it looks more like proportions ?

[Thomas Rodriguez]

Yes, the results should be proportions of a population becuase one of the theorems for the SEIR model is N (Our Total Population) = S + E + I + R = 1. Hence the results from each repsective population should be proportions of 1 to fulfill the theorem shown.

Question:

Would [ t , x ] = ode23( odefun, tspan, x0) contain a remoldeled Runge Kutta method to solve for the next vectors of x?

[Walter Roberson]

I do not understand what you mean about "remolded" RK method? You have differential equations, and the ode*() functions would permit you to do numeric analysis over time.

[Thomas Rodriguez]

Sorry I misinterpreted your comment. But I reorganinzed the code to the following:

clc;clear;close all;
r0 = 3;
SSm = r0;%Social Measures
Tinc = 0.5; %Time of Incubation ~ 2 days
Tinf = 0.1; %Time of Infected ~ 10 days
t0 = 0; %Initial Time - Given
tf = 300;%Final Time - Given
N = 300; %Number of Points (or "Day Intervals") % M = N 
tol = 10^(-6); %Tolerance

h = (tf-t0)/N;%Step Size
t = [t0:h:tf]';%Discretized time step

S0 = 0.8; %Initial amount of Susceptibles
E0 = 0.2; %Initial amount of Incubants
I0 = 0.0; %Initial amount of Infectants
R0 = 0.0; %Initial amount of Recovered

%Initial Values X(0) - Matrix (4x1)
x0 = [S0;E0;I0;R0];
x = zeros(4,N); x(:,1) = x0 ; 
Sm = zeros(1,N); Sm(1) = SSm;
% f(t,x) = dx/dt
%f = @(t,x) [(-R*x(3)*x(1)/Tinf);(R*x(3)*x(1)/Tinf)-(x(2)/Tinc);(x(2)/Tinc)-(x(3)/Tinf);(x(3)/Tinf)];
dS = @(s,e,i,r) (-SSm.*i.*s/Tinf);
dE = @(s,e,i,r) (SSm.*i.*s/Tinf)-(e/Tinc);
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);
dR = @(s,e,i,r) (i/Tinf);
ODE = @(t,x) [dS(x(1),x(2),x(3),x(4));
    dE(x(1),x(2),x(3),x(4));
    dI(x(1),x(2),x(3),x(4));
    dR(x(1),x(2),x(3),x(4));];

[t,x] = OdeSolver(ODE,t,x);


%Tables:
fprintf('\n')
fprintf('----------------------------------------------------------------------------------------------------------------------\n')
fprintf('----------------------------------------------------------------------------------------------------------------------\n')
fprintf('                                SEIR Model Stats                                     \n')
fprintf('----------------------------------------------------------------------------------------------------------------------\n')
fprintf('iter(k)  Days         Susceptible       Incubants         Infectants        Recovered\n')
for i = 1:N
    fprintf('%2.0f   %7.6e    %13.6e    %13.6e     %13.6e     %13.6e\n', i, t(i), x(1,i), x(2,i), x(3,i), x(4,i))
    fprintf('----------------------------------------------------------------------------------------------------------------------\n')
end
S = x(1,:);
E = x(2,:);
I = x(3,:);
R = x(4,:);
%Plotting
figure(1);hold on; box on; grid on;
plot(t,S,t,E,t,I,t,R)
xlabel('Time (Days)')
ylabel('dx/dt')
title('SEIR Model')
legend( 'Susceptible', 'Incubants', 'Infected', 'Recovered')
hold off;
%% 
function [t,x] = OdeSolver(ODE,t,x)
%Given Initial Values
r0 = 3;
SSm = r0;%Social Measures
Tinc = 0.5; %Time of Incubation ~ 2 days
Tinf = 0.1; %Time of Infected ~ 10 days
t0 = 0; %Initial Time - Given
tf = 300;%Final Time - Given
N = 300; %Number of Points (or "Day Intervals") % M = N 
h = (tf-t0)/N;%Step Size
%ODEs for SEIR
dS = @(s,e,i,r) (-SSm.*i.*s/Tinf);%dS/dt - Susceptible Rate
dE = @(s,e,i,r) (SSm.*i.*s/Tinf)-(e/Tinc);%dE/dt - Incubant Rate
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);%dI/dt - Infectants Rate
dR = @(s,e,i,r) (i/Tinf);%dR/dt - Recovered Rate
for i = 1:N
    %k1 = h*f(tn,yn)
    s1 = h*dS( x(1,i), x(2,i), x(3,i), x(4,i));
    e1 = h*dE( x(1,i), x(2,i), x(3,i), x(4,i));
    i1 = h*dI( x(1,i), x(2,i), x(3,i), x(4,i));
    r1 = h*dR( x(1,i), x(2,i), x(3,i), x(4,i));
    %k2 = h*f(tn + (h/2),yn + k1.*(h/2))
    s2 = h*dS( x(1,i) + 0.5*s1, x(2,i) + 0.5*e1, x(3,i) + 0.5*i1, x(4,i) + 0.5*r1);
    e2 = h*dE( x(1,i) + 0.5*s1, x(2,i) + 0.5*e1, x(3,i) + 0.5*i1, x(4,i) + 0.5*r1);
    i2 = h*dI( x(1,i) + 0.5*s1, x(2,i) + 0.5*e1, x(3,i) + 0.5*i1, x(4,i) + 0.5*r1);
    r2 = h*dR( x(1,i) + 0.5*s1, x(2,i) + 0.5*e1, x(3,i) + 0.5*i1, x(4,i) + 0.5*r1);
    %k3 = h*f(tn + (h/2),yn + k2.*(h/2))
    s3 = h*dS( x(1,i) + 0.5*s2, x(2,i) + 0.5*e2, x(3,i) + 0.5*i2, x(4,i) + 0.5*r2);
    e3 = h*dE( x(1,i) + 0.5*s2, x(2,i) + 0.5*e2, x(3,i) + 0.5*i2, x(4,i) + 0.5*r2);
    i3 = h*dI( x(1,i) + 0.5*s2, x(2,i) + 0.5*e2, x(3,i) + 0.5*i2, x(4,i) + 0.5*r2);
    r3 = h*dR( x(1,i) + 0.5*s2, x(2,i) + 0.5*e2, x(3,i) + 0.5*i2, x(4,i) + 0.5*r2);
    %k4 = h*f(tn + h,yn + h.*k3)
    s4 = h*dS( x(1,i) + s3, x(2,i) + e3, x(3,i) + i3, x(4,i) + r3);
    e4 = h*dE( x(1,i) + s3, x(2,i) + e3, x(3,i) + i3, x(4,i) + r3);
    i4 = h*dI( x(1,i) + s3, x(2,i) + e3, x(3,i) + i3, x(4,i) + r3);
    r4 = h*dR( x(1,i) + s3, x(2,i) + e3, x(3,i) + i3, x(4,i) + r3);
    %yn+1 = yn + (1/6)*h*(k1+2*k2+2*k3+k4)
    x(1,i+1) = x(1,i) + (s1 + 2*s2 + 2*s3 + s4)/6;
    x(2,i+1) = x(2,i) + (e1 + 2*e2 + 2*e3 + e4)/6;
    x(3,i+1) = x(3,i) + (i1 + 2*i2 + 2*i3 + i4)/6;
    x(4,i+1) = x(4,i) + (r1 + 2*r2 + 2*r3 + r4)/6;
end
end

After running the code I got the following results:

----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------------
                                SEIR Model Stats                                     
----------------------------------------------------------------------------------------------------------------------
iter(k)  Days         Susceptible       Incubants         Infectants        Recovered
 1   0.000000e+00     8.000000e-01     2.000000e-01      0.000000e+00      0.000000e+00
----------------------------------------------------------------------------------------------------------------------
 2   1.000000e+00     3.570400e+03    -3.585533e+03     -6.533333e+00      2.266667e+01
----------------------------------------------------------------------------------------------------------------------
 3   2.000000e+00     3.548099e+28    -3.548099e+28     -6.267583e+16     -6.290141e+10
----------------------------------------------------------------------------------------------------------------------
 4   3.000000e+00    3.153960e+201    -3.153960e+201     -5.921899e+123     -5.917697e+76
----------------------------------------------------------------------------------------------------------------------
 5   4.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------
 6   5.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------
 7   6.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------
 8   7.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------
 9   8.000000e+00              NaN              NaN               NaN               NaN
 

From the results, the code works fine until the 5th iteration. You mentioned that it could be the dE or (both) dI ODEs. How could it be these equations that create this effect?

[Walter Roberson]

dS = @(s,e,i,r) (-SSm.*i.*s/Tinf);
dE = @(s,e,i,r) (SSm.*i.*s/Tinf)-(e/Tinc);
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);
dR = @(s,e,i,r) (i/Tinf);
ODE = @(t,x) [dS(x(1),x(2),x(3),x(4));
    dE(x(1),x(2),x(3),x(4));
    dI(x(1),x(2),x(3),x(4));
    dR(x(1),x(2),x(3),x(4));];

It does not make sense to combine that with

%ODEs for SEIR
dS = @(s,e,i,r) (-SSm.*i.*s/Tinf);%dS/dt - Susceptible Rate
dE = @(s,e,i,r) (SSm.*i.*s/Tinf)-(e/Tinc);%dE/dt - Incubant Rate
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);%dI/dt - Infectants Rate
dR = @(s,e,i,r) (i/Tinf);%dR/dt - Recovered Rate
for i = 1:N

There is no difference between the anonymous functions, so why would you do the assignments in both locations?

If you need the parts broken out individually, then in the outer routine define

  ODE = {dS, dE, dI, dR}

and inside the function,

  dS = ODE{1}; dE = ODE{2}; dI = ODE{3}; dR = ODE{4};

Anyhow, earlier, I used your dS dE dI dR the way you had provided them, and I ran your code, and I saw the same kind of NaN that you describe. I then defined a function like the ODE = @(t,x) etc that you define here, and used it together with Mathwork's ode23() call. ode23() did not encounter nans when passed the function. I could be mistaken when I thought your code looked like an implementation of RK23 (for example perhaps your implemation was more like ode15()). It looked to me as if the fact that your implementation of what appears to be RK23 encounters nans but ode23() does not for the same functions, suggests to me that your RK implementation is wrong. I suspect that if you test with simple functions that you will not see correct results from your RK implementation.

The fact that I saw negative and positive for the second and third functions suggests to me that your equations are also wrong -- but that the equations being wrong does not explain why you got the nan.

Answer 2

[the cyclist]

In this line

S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6;

the variables s1, s2, s3, and s4 are all row vectors with 300 elements.

The left-hand side is a single element.

I haven't tried to figure out what you intended to do, but that's the basic problem.

[Thomas Rodriguez]

My intentions for the algorthim is to have it calculate s1,s2,s3, and s4. Then to compute the next iteration of the loop it would use the S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6 and storing those "S" values into the S array for later access.

[Thomas Rodriguez]

%Initial Values X(0) - Matrix (4x1)
%S = zeros(1,N);
S(1) = S0;
%E = zeros(1,N); 
E(1) = E0;
%I = zeros(1,N);
I(1) = I0;
%R = zeros(1,N); 
R(1) = R0;

So I've adjusted my code a bit to instead of calling upon a full array of zeros that is length N. I just treated it as an empty array. Results seem to appear but not my desire ones:

----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------------
                                SEIR Model Stats                                     
----------------------------------------------------------------------------------------------------------------------
iter(k)  Days         Susceptible       Incubants         Infectants        Recovered
 1   0.000000e+00     8.000000e-01     2.000000e-01      0.000000e+00      0.000000e+00
----------------------------------------------------------------------------------------------------------------------
 2   1.000000e+00     8.000000e-01     6.666667e-02     -1.453333e+01      1.466667e+01
----------------------------------------------------------------------------------------------------------------------
 3   2.000000e+00     8.000000e-01     2.222222e-02     -4.234044e+03      4.234222e+03
----------------------------------------------------------------------------------------------------------------------
 4   3.000000e+00     8.000000e-01     7.407407e-03     -1.232109e+06      1.232109e+06
----------------------------------------------------------------------------------------------------------------------
 5   4.000000e+00     8.000000e-01     2.469136e-03     -3.585436e+08      3.585436e+08
----------------------------------------------------------------------------------------------------------------------
 6   5.000000e+00     8.000000e-01     8.230453e-04     -1.043362e+11      1.043362e+11
----------------------------------------------------------------------------------------------------------------------
 7   6.000000e+00     8.000000e-01     2.743484e-04     -3.036183e+13      3.036183e+13
----------------------------------------------------------------------------------------------------------------------

As you can see the Susceptible (aka 'S') is not creating new values for the next iterations, but the other populations E,I,and R are. Could this be from the layout of the loop for S more specifically the statement: S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6?

[the cyclist]

I notice that in these anonymous functions, such as

dS = @(s,e,i,r) (-R.*i.*s/Tinf);
dE = @(s,e,i,r) (R.*i.*s/Tinf)-(e/Tinc);

the input argument is r, but the expression evaluated has R.

I'm guessing that that was not intentional. It looks like dS will therefore always be zero, instead of updating with the other expressions.

[Thomas Rodriguez]

R is a scalar values defined at the beginning

r0 = 3;
R = r0;%For Now

The purpose of R is to add propoagation in the ODE by entering as a scalar value but for future calculations it will be defined as a function component of time. But for now it is a scalar value. Also the reason that 'r' is not present in dS, its because the amount of "Recovered"(R) individuals does not effect the amount of Susceptible(S).

[Walter Roberson]

Then to compute the next iteration of the loop it would use the S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6

Use cell indexing instead of () indexing.

S{i+1} = S{i} + (s1 + 2*s2 + 2*s3 + s4)/6

[Thomas Rodriguez]

From using

S{i+1} = S{i} + (s1 + 2*s2 + 2*s3 + s4)/6;

I recieved the following errors:

"Brace indexing is not supported for variables of this type.

Error in SEIR_RK4 (line 59)

S{i+1} = S{i} + (s1 + 2*s2 + 2*s3 + s4)/6;"

Then I changed the array S into a cell array:

S{1} = S0;

Then I got the following results:

----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------------
                                SEIR Model Stats                                     
----------------------------------------------------------------------------------------------------------------------
iter(k)  Days         Susceptible       Incubants         Infectants        Recovered
 1   0.000000e+00     8.000000e-01     2.000000e-01      0.000000e+00      0.000000e+00
----------------------------------------------------------------------------------------------------------------------
 2   1.000000e+00     3.570400e+03    -3.585533e+03     -6.533333e+00      2.266667e+01
----------------------------------------------------------------------------------------------------------------------
 3   2.000000e+00     3.548099e+28    -3.548099e+28     -6.267583e+16     -6.290141e+10
----------------------------------------------------------------------------------------------------------------------
 4   3.000000e+00    3.153960e+201    -3.153960e+201     -5.921899e+123     -5.917697e+76
----------------------------------------------------------------------------------------------------------------------
 5   4.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------
 6   5.000000e+00              NaN              NaN               NaN               NaN
----------------------------------------------------------------------------------------------------------------------

The algorthim only did successfully 4 iterations and the rest of the 300 exploded to infinity it seems.

I also got the following error from attempting to plot the data:

"Error using plot

Invalid data argument.

Error in SEIR_RK4 (line 79)

plot(t,S,t,E,t,I,t,R)"