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);

Categories

Find more on Modeling in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!