Azzera filtri
Azzera filtri

This code runs except the graph and results are incredibly wrong

1 visualizzazione (ultimi 30 giorni)
This code is running except the graph and the results are really wrong. I will provide the correct graph as a refference. I have checked all the equations and parameters but I cant seem to find why the answers are going to values that are *10^26
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; Bi];
[t,y] = rk4(@PSLIR, tspan, y0, Bmax, uL, uI, e, T, tspan, Ap);
figure
hold on
plot(t,y(1,:),'k')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
plot(t,y(4,:),'y')
plot(t,y(5,:),'m')
legend("S","L","I","R","P")
%% 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(5);
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 = odeFunc(t(n), y(:, n), Bmax, uL, uI, e, T, day, Ap);
k2 = odeFunc(t(n) + 0.5*h, y(:,n) + 0.5.*k1*h, Bmax, uL, uI, e, T, day, Ap);
k3 = odeFunc(t(n) + 0.5*h, y(:,n) + 0.5*k2*h, Bmax, uL, uI, e, T, day, Ap);
k4 = odeFunc(t(n) + h, y(:,n) + k3*h, Bmax, uL, uI, e, T, day, Ap);
y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))/6;
end
end
  1 Commento
CONNOR
CONNOR il 1 Dic 2023
I apoligise I hadnt linked the files. Now the info file and the correct graph file have been linked

Accedi per commentare.

Risposta accettata

David Goodmanson
David Goodmanson il 3 Dic 2023
Modificato: David Goodmanson il 3 Dic 2023
Hi Connor,
The reason your code is going to 10^26 is that you are missing an important factor on the last line of rk4, which is the width of the interval. You should have
y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))*h/6;
(running a test case on rk4 before trying to use it in the larger code might well have saved you time overall). After that change, you get a finite solution which takes care of the 10^26 effect.
The solution does not agree with the plot, and to agree with that plot you would have to adjust some of the starting values. Also it appears that the code builds up y0 as P,S,L,I,B but plots out y as S,L,I,R,P

Più risposte (0)

Categorie

Scopri di più su Startup and Shutdown in Help Center e File Exchange

Prodotti


Release

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by