I don't know why this code is saying, "Arrays have incompatible sizes for this operation."

1 view (last 30 days)
This code has been giving me a lot of trouble trying to get it to work and now it is saying that the arrays are incompatable but I dont know why that is and where to look to fix that.
close all, clear, clc
load('EnvironmentalForcing.mat')
Bmax = 1;
uL_min = 1/6;
uI = 1/10;
e = 0.001;
Ap = 5000;
Pi = 930.27249;
Si = Pi/Ap;
Li = 0.01*Si;
Ii = 0;
Ri = uI*Ii;
Bi = 1e-4;
T_beta = zeros(1, length(tspan));
uL = zeros(1, length(tspan));
for i=1:length(tspan)
if (T(i) <= 0) || (T(i)>=35)
T_beta(i) = 0;
else
T_beta(i) = 0.00241*T(i)^2.06737 * (35-T(i))^0.72859;
end
end
for i=1:length(tspan)
uL(i) = 1/sum(T_beta(1:i));
j = 0;
while 1/uL(i) > 1/uL_min
j = j+1;
uL(i) = 1/sum(T_beta(1+j: i));
end
end
% Te = -0.35068 + 0.10789.*T -0.00214.*T.^2;
% for i = 1:length(T)
% if T(i)>0 && T(i)<35
% Tb(i) = (0.000241*(T(i)^2.06737))*((35-T(i))^0.72859);
% end
% end
% j = 1;
% for i=1:length(t)
% uL(i) = sum(Tb(j:i));
% while uL(i) > uL_min
% j = j+1;
% uL(i) = sum(Tb(j:i));
% end
% end
% B = Bmax*Tb;
% p = [Bmax, uL, uI, e, Te];
y0 = [Pi; Si; Li; Ii; Ri; Bi];
[t,y] = rk4(@PSLIR, tspan, y0, Bmax, uL, uI, e, T, tspan, Ap);
%% Functions
function [dydt] = PSLIR(t_index, y0, Bmax, uL, uI, e, T, day, Ap)
if (isa(t_index, 'int8')) &&(t_index >= 2)
t_index_rounded = round(t_index);
elseif (t_index <2)
t_index_rounded = 1;
else
t_index_rounded = round(t_index)-1;
end
if (t_index_rounded > 0) && (t_index_rounded <= length(T))
Te = -0.35068 + 0.10789.*T(t_index_rounded) -0.00214.*T(t_index_rounded).^2;
if (T(t_index_rounded) <=0) || (T(t_index_rounded) >=35)
T_beta = 0;
else
T_beta = (0.000241*(T(t_index_rounded)^2.06737))*((35-T(t_index_rounded))^0.72859);
end
B = Bmax.*T_beta;
else
error('Cant compute')
end
%assign variables
P = y0(1);
S = y0(2);
L = y0(3);
I = y0(4);
R = y0(5);
Pb = y0(6);
dPbdt = (0.1724.*Pb-0.0000212.*Pb^2).*Te;
dPdt = ((1.33.*day(t_index_rounded)).*Te)+dPbdt;
dSdt = (-B.*S.*I)+dPdt./Ap;
dLdt = (B.*S.*I)-((uL(t_index_rounded).^-1).*L)+e;
dIdt = ((uL(t_index_rounded).^-1).*L)-((uI.^-1).*I);
dRdt = (uI.^-1).*I;
dydt = [dPdt; dSdt; dLdt; dIdt; dRdt];
end
function [t, y] = rk4(odeFunc, tspan, y0, Bmax, uL, uI, e, T, day, Ap)
% N = length(tspan);
% q = length(y0);
% t0 = tspan(2);
% h = tspan(2) - tspan(1);
% t = zeros(N, 1);
% y = zeros(q, N);
% t(1) = t0;
% y(:, 1) = y0;
N = length(tspan);
t = tspan;
y = zeros(length(y0),N);
y(:,1) = y0(:);
for n=1:N-1
h = tspan(n+1)-tspan(n);
k1 = h*odeFunc(n, y(:, n), Bmax, uL, uI, e, T, day, Ap);
k2 = h*odeFunc(n + 0.5, y(:,n) + 0.5.*k1, Bmax, uL, uI, e, T, day, Ap);
k3 = h*odeFunc(n + 0.5, y(:,n) + 0.5*k2, Bmax, uL, uI, e, T, day, Ap);
k4 = h*odeFunc(n + h, y(:,n) + k3, Bmax, uL, uI, e, T, day, Ap);
y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))/6;
end
end

Accepted Answer

Torsten
Torsten on 1 Dec 2023
Edited: Torsten on 1 Dec 2023
Your vector of initial conditions is 6x1
y0 = [Pi; Si; Li; Ii; Ri; Bi];
but you only return a 5x1 vector in PSLIR
dydt = [dPdt; dSdt; dLdt; dIdt; dRdt];
Further, your first arguments in the calls to PSLIR are wrong (at least according to Runge-Kutta 4th order).
They must be t(n), t(n)+0.5*h, t(n)+0.5*h, t(n)+h in the calculation of k1,k2,k3 and k4 instead of n, n+0.5, n+0.5, n+h.

More Answers (0)

Categories

Find more on Dates and Time in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!