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

1 次查看(过去 30 天)
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

采纳的回答

Torsten
Torsten 2023-12-1
编辑:Torsten 2023-12-1
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.

更多回答(0 个)

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by