This code runs except the graph and results are incredibly wrong

1 次查看（过去 30 天）
CONNOR
CONNOR 2023-12-1
编辑： David Goodmanson 2023-12-3
Ran in:
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

请先登录，再回答此问题。

采纳的回答

David Goodmanson
David Goodmanson 2023-12-3
编辑：David Goodmanson 2023-12-3
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Startup and Shutdown 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by