Index exceeds array bounds
显示 更早的评论
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
采纳的回答
Walter Roberson
2020-8-1
0 个投票
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.
3 个评论
Thomas Rodriguez
2020-8-1
Walter Roberson
2020-8-1
x(1,1) = S0;
x(2,1) = E0;
x(3,1) = I0;
x(4,1) = R0;
x is a column vector of length 4
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)];
f expects a two parameters to be passed in, with the second one named x for the purpose of processing f. Inside f, x refers to whatever was passed as the second parameter, not to the column vector of length 4. The function f expects the second argument to be an array with at least 3 elements (but does not need 4 elements for it.)
[G, Grad] = Gn(x,t,h,f,R,Tinc,Tinf);
Gn is passed the complete column vector of length 4, and also is passed the anonymous function f.
function [G, Grad] = Gn(x,t,h,f,R,Tinc,Tinf)
G = @(x) x - x - (h/2)*( f(t,x) + f(t+h,(x+h)) ); %Gn(x)
G is defined as an anonymous function that takes one parameter. The parameter is named x for the purpose of processing G. Inside G, x refers to whatever was passed as the parameter, not to whatever was passed as the first parameter of Gn, and not to the column vector of length 4. For reasons that are not clear, it subtracts the parameter from itself, which gives an array of zeros the same size as the parameter. That appears to be a waste of time unless the purpose is to provide shaping information for how to replicate the rest of the expression if the rest of the expression produces something of different size. During the rest of the expression, the parameter x is passed to f. Recall that f expects its second parameter to be at least 3 elements, so G will end up failing if the parameter that is passed to it is not at least length 3.
%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);
The scalar value x(i) and the anonymous functions G and Grad are passed to Newtons_Method
function [x] = Newtons_Method(x,G,Grad)
which happens to receive them under the same name. But remember that the x received here will be the scalar known as x(i) in the calling function, not a vector.
while( (abs(norm(G(x))) > tol && loop_0 <= MaxLoop && norm(Grad(x)) > 0) )
The scalar value known locally as x will be passed to G.
G = @(x) x - x - (h/2)*( f(t,x) + f(t+h,(x+h)) ); %Gn(x)
G will pass it to f
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)];
f will try to access index 3 of the scalar, which will fail.
%Newtons Method: X_n+1 = X_n - inv(Grad(X_n))*G(X_n)
That notation means that the (n+1)'s "generation" of the vector X is determined by applying that series of operations to the n'th "generation" of the vector X . The notation does not mean that the n+1'st component of vector X is to be calculated that way. You have to read it like
X{n+1} = X{n} - inv(Grad(X{n})) * G(X{n})
Thomas Rodriguez
2020-8-1
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
