MATLAB: Index exceeds array bounds

indexingmatrix arraynonlinear

Hello all! I'm developing an SEIR Model code that is solved through the Trapeziodal Method Implicitly but uses Newton's method to find the next vector. I keeping getting the error "Index exceeds array bounds:
Error in
SEIR_Model>@(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)]
Error in SEIR_Model>@(x)x-x-(h/2)*(f(t,x)+f(t+h,(x+h))) (line 75)
G = @(x) x – x – (h/2)*( f(t,x) + f(t+h,(x+h)) ); %Gn(x)
Error in SEIR_Model>Newtons_Method (line 101)
while( (abs(norm(G(x))) > tol && loop_0 <= MaxLoop && norm(Grad(x)) >
0) )
Error in SEIR_Model (line 38)
x(:,(i+1)) = Newtons_Method(x(i),G,Grad);"
Code Shown Below:
clc;clear;close all;
r0 = 3;
R = r0;%For Now
Tinc = 0.1; %Time of Incubation
Tinf = 0.2; %Time of Infected
t0 = 0; %Initial Time - Given
tf = 300;%Final Time - Given
N = 299; %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)
x = zeros(4,N);%Creating matrix 4xM
x0 = [S0;E0;I0;R0];
x(1,1) = S0;
x(2,1) = E0;
x(3,1) = I0;
x(4,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)];
 %Trapezoidal Method: X_n+1 = X_n + (h./2)*( f(t_n+1,X_n+1) + f(t_n,X_n) )
%[G, Grad] = Gn(x,t,h,f,R,Tinc,Tinf);
for i = 1:N
    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)];
    [G, Grad] = Gn(x,t,h,f,R,Tinc,Tinf);
    %Newtons Method: X_n+1 = X_n - inv(Grad(X_n))*G(X_n)
    %Adding the new values into the 4xM matrix 'x'
    x(:,(i+1)) = Newtons_Method(x(i),G,Grad);
end
%Extracting Values of each Population (S,E,I,R) by calling the values for
%each of the Rows through the calling below
S = []; S = [S x(1,:)];
E = []; E = [E x(2,:)];
I = []; I = [I x(3,:)];
R = []; R = [R x(4,:)];
%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), S(i), E(i), I(i), R(i))
    fprintf('----------------------------------------------------------------------------------------------------------------------\n')
end
%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 [G, Grad] = Gn(x,t,h,f,R,Tinc,Tinf)
% X_n+1 = X_n + (h./2)*( f(t_n+1,X_n+1) + f(t_n,X_n) )
% Creating Implicit form of the Trapeziodal Method
% Gn(X) = x - x - (1/2)*(f(t_n+1,X_n+1) + f(t_n,X_n)) 
G = @(x) x - x - (h/2)*( f(t,x) + f(t+h,(x+h)) ); %Gn(x)
G0 = eye(4);
Grad = @(x) G0-(h/2)*[(-R*x(3,end))/(Tinf) 0 (-R*x(1,end))/(Tinf) 0;
    (R*x(3,end)/Tinf)-(x(2,end)/Tinc) (R*x(3,end)*x(1,end)/Tinf)-(1/Tinc) (R*x(1,end)/Tinf)-(x(2,end)/Tinc) 0;
    0 (1/Tinc) -(1/Tinf) 0;
    0 0 (1/Tinf) 0]; %Gradient of Gn(x)- 4x4 matrix composed of the following
% [Equation 1: dS dE dI dR; 
% Equation 2:dS dE dI dR; 
% Equation 3:dS dE dI dR; 
% Equation 4:dS dE dI dR;]
%Printing Table of Gn(X)
%fprintf('Gn(x) - Trapezoidal Method Reformulated into an Implicit Form:\n')
%disp(G)
%fprintf('Gradient (Jacobian) of Gn(x) Form:\n')
%disp(Grad)
end
%% 
function [x] = Newtons_Method(x,G,Grad)
%Calling upon the Gn(x) & Gradient function
tol = 10^(-6); %Tolerance
loop_0 = 1; %Initial Loop count
MaxLoop = 300; %Maximum Number of Iterations
%&& norm(G(x)) > 0
while( (abs(norm(G(x))) > tol && loop_0 <= MaxLoop && norm(Grad(x)) > 0) )
    gx = G(x);
    Gdx = Grad(x);
    x = x - gx/Gdx;
    loop_0 = loop_0 + 1;
end
%Discussion of Algorthim: 
% Tolerance (tol) was chosen by user which dictates how minimal the error of our Newton's Method with the tolernace
% nGrad > 0 - This was created as a condition to make sure that their is smoothness for the Gradient or dG/dx.
% The nG & nGrad are calculated because it converts the matrices of both the Gradient & Gn(x) into a scalar values by taking the norm
end

Best Answer

Your f requires that x is a vector.
You call Newtons_Method passing in a scalar for the first value. The function receives it as x. The function calls G with that scalar x. G passes the scalar x to f. But f needs vector x.
Related Question