I apoligise I hadnt linked the files. Now the info file and the correct graph file have been linked

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

1 view (last 30 days)

Show older comments

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

### Accepted Answer

David Goodmanson
on 3 Dec 2023

Edited: David Goodmanson
on 3 Dec 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

##### 0 Comments

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!