i need to find error in code because I it's increasing but it should be decresing

1 view (last 30 days)
i need to find error in code because I it's increasing but it should be decresing
the equation of Sir model in runge kutta (hune's method)
this is laws
i get the constants on this picture and laws
this my code
h=1.5;
a=4*10^-6; %a=alfa
y=0.017; %gamaa
S=zeros(22,1);
S(1)=54000;
I=zeros(22,1);
I(1)=757;
R=zeros(22,1);
R(1)=0;
fprintf('\n')
disp(' heun`s method ')
disp('-----------------------------------------------------------------------------------------------')
disp(' t S I R')
disp('-----------------------------------------------------------------------------------------------')
fprintf('\n')
for i=1:22
k1=-a*S(i)*I(i);
l1=a*S(i)*I(i)-y*I(i);
m1=y*I(i);
k2=-a*(S(i)+h*k1)*(I(i)+h);
l2=a*(S(i)+h*k1)*(I(i)+h)-y*(I(i)+h);
m2=y*(I(i)+h);
S(i+1)=S(i)+(h/2)*(k1+k2);
I(i+1)=I(i)+(h/2)*(l1+l2);
R(i+1)=R(i)+(h/2)*(m1+m2);
fprintf('%10f %17f %17f %17f\n',i,S(i+1),I(i+1),R(i+1))
end
t=[1:22];
hold on
plot(S(t),t,'-b')
plot(I(t),t,'--r')
plot(R(t),t,'-g')

Answers (2)

Mathieu NOE
Mathieu NOE on 28 Apr 2022
hi
seems you have swapped the x and y data when you do the plot
after the fix we get this :
full code a bit improved :
h=1.5;
alfa=4*10^-6; %a=alfa
gamma=0.017; %gamma
n = 22; % time steps (samples)
S=zeros(n,1);
I=zeros(n,1);
R=zeros(n,1);
S(1)=54000;
I(1)=757;
R(1)=0;
fprintf('\n')
disp(' heun`s method ')
disp('-----------------------------------------------------------------------------------------------')
disp(' t S I R')
disp('-----------------------------------------------------------------------------------------------')
fprintf('\n')
for i=1:n
k1=-alfa*S(i)*I(i);
l1=alfa*S(i)*I(i)-gamma*I(i);
m1=gamma*I(i);
k2=-alfa*(S(i)+h*k1)*(I(i)+h);
l2=alfa*(S(i)+h*k1)*(I(i)+h)-gamma*(I(i)+h);
m2=gamma*(I(i)+h);
S(i+1)=S(i)+(h/2)*(k1+k2);
I(i+1)=I(i)+(h/2)*(l1+l2);
R(i+1)=R(i)+(h/2)*(m1+m2);
fprintf('%10f %17f %17f %17f\n',i,S(i+1),I(i+1),R(i+1))
end
t=[1:n];
plot(t,S(t),'-b',t,I(t),'--r',t,R(t),'-g')
legend('S','I','R');

James Tursa
James Tursa on 3 May 2022
Edited: James Tursa on 3 May 2022
You need to use all the derivatives at the first point to propagate for the initial guess at the second point. I.e., these lines
k2=-a*(S(i)+h*k1)*(I(i)+h);
l2=a*(S(i)+h*k1)*(I(i)+h)-y*(I(i)+h);
m2=y*(I(i)+h);
should be this instead, where (I(i)+h) is changed to (I(i)+h*l1):
k2=-a*(S(i)+h*k1)*(I(i)+h*l1);
l2=a*(S(i)+h*k1)*(I(i)+h*l1)-y*(I(i)+h*l1);
m2=y*(I(i)+h*l1);

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!