# This code runs except the graph and results are incredibly wrong

조회 수: 2 (최근 30일)
CONNOR 2023년 12월 1일
편집: David Goodmanson 2023년 12월 3일
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
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
CONNOR 2023년 12월 1일
I apoligise I hadnt linked the files. Now the info file and the correct graph file have been linked

댓글을 달려면 로그인하십시오.

### 채택된 답변

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

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile 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