Index exceeds the number of array elements. Index must not exceed 3. HELP

model code:

clear all;

close all;

S = 99;

I = 1;

R = 0;

N = 100; %Total population

beta= 0.1; % birth rate

alpha= 0.1; % infection person to person rate

lambda= 0.3; % infection by water rate

vac= 0.05; % recovery by vaccination rate

d= 0.03; % death rate

gamma= 0.8; % recovery rate

c= 0.9; % rate of contamination

m= 0.4; % rate of decay of V. cholerae

B= 0.0; % initial concentration of V. cholerae

t_f = 500; %Ending time of simulation

Q = [beta alpha lambda vac d gamma c m B];

[t,y]=ode45('cholera_de',[0:t_f/100:t_f],[S I R]',[],Q);

figure(1)

plot(t,y(:,1),'k-',t,y(:,2),'r--',t,y(:,3),'b:');

xlabel('\bf Time (days)');

ylabel('\bf Number of People by Category');

legend('S','I','R');

z=y(end,:)'

SN=y(:,1)/N;

IN=y(:,2)/N;

figure(2)

plot(IN,SN);

xlabel('\bf I/N');

ylabel('\bf S/N');

r0=beta/beta

[maxIN,y_maxtime]=max(y(:,2)/N);

maxIN

maxtime=y(y_maxtime)

eqIN=y(100,2)/N;

eqIN

function code:

function dy=cholera_de(t,Y,flag,Q)

beta= Q(1);

alpha= Q(2);

lambda= Q(3);

vac= Q(4);

d= Q(5);

gamma= Q(6);

c= Q(7);

m= Q(8);

S= Y(1);

I= Y(2);

R= Y(3);

B=Y(4);

N= S+I+R;

dy(1,1)= beta - alpha*I - lambda*B - vac*S - d;

dy(2,1)= alpha*I + lambda*B - d - gamma*I;

dy(3,1)= gamma*I + vac*S - d;

dy(4,1)= c*I - m*B;

Gives me this error:

Index exceeds the number of array elements. Index must not exceed 3.

Error in cholera_de (line 15)

B=Y(4);

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in cholera_model (line 23)

[t,y]=ode45('cholera_de',[0:t_f/100:t_f],[S I R]',[],Q);

[Chuguang Pan]

It seems that your initial condition is [S I R]', which has three elements. However, the cholera_de function need four states Y(1), Y(2), Y(3), Y(4). You should keep number of initial states the same as states transfered to cholera_de function.

Answer 1

[Sam Chak]

Hi @Elizabeth Pickford

Few places require fixes. Check out the annotated lines.

S0      = 99;
I0      = 1;
R0      = 0;
N0      = S0 + I0 + R0;     % Total population              % <-- fix it here

beta    = 0.1;              % birth rate
alpha   = 0.1;              % infection person to person rate
lambda  = 0.3;              % infection by water rate
vac     = 0.05;             % recovery by vaccination rate
d       = 0.03;             % death rate
gamma   = 0.8;              % recovery rate
c       = 0.9;              % rate of contamination
m       = 0.4;              % rate of decay of V. cholerae
B       = 0.0;              % initial concentration of V. cholerae
Q       = [beta alpha lambda vac d gamma c m B];

t_f     = 500;              % Ending time of simulation
tspan   = [0:t_f/100:t_f];                                  % <-- fix it here
y0      = [S0; I0; R0; N0];                                 % <-- fix it here
[t, y]  = ode45(@(t, y) cholera_de(t, y, Q), tspan, y0);    % <-- fix it here

figure(1)
plot(t,y(:,1),'k-',t,y(:,2),'r--',t,y(:,3),'b:'); grid on
xlabel('\bf Time (days)');
ylabel('\bf Number of People by Category');
legend('S', 'I', 'R');

z   = y(end,:)'

z = 4x1

14.2061 -0.9011 91.6949 -2.0108

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

N   = y(:,1) + y(:,2) + y(:,3);                             % <-- fix it here
SN  = y(:,1)./N;                                            % <-- fix it here
IN  = y(:,2)./N;                                            % <-- fix it here

figure(2)
plot(IN, SN); grid on
xlabel('\bf I/N');
ylabel('\bf S/N');

r0      = beta/beta;
[maxIN, y_maxtime] = max(y(:,2)/N);
maxIN;
maxtime = y(y_maxtime);
eqIN    = y(100,2)/N;
eqIN;

function dy = cholera_de(t, y, Q)                           % <-- fix it here
    %% parameters
    beta    = Q(1);
    alpha   = Q(2);
    lambda  = Q(3);
    vac     = Q(4);
    d       = Q(5);
    gamma   = Q(6);
    c       = Q(7);
    m       = Q(8);
    
    %% definitions
    S       = y(1);
    I       = y(2);
    R       = y(3);
    B       = y(4);
    N       = S + I + R;
    
    %% differential equations
    dy(1,1) = beta    - alpha*I  - lambda*B - vac*S - d;
    dy(2,1) = alpha*I + lambda*B - gamma*I  - d;
    dy(3,1) = gamma*I + vac*S    - d;
    dy(4,1) =     c*I - m*B;
end