1 view (last 30 days)
Surath Ghosh on 3 Dec 2021
Commented: Surath Ghosh on 3 Dec 2021
I have written the code but it is unable to produce correct result. Can anybody help me to find the errors in the above code?
Walter Roberson on 3 Dec 2021
No obvious error in your output? I do not see any error message ?
As outside observers who have not been given documentation, we must assume that this is what your code was designed to calculate.
Otherwise we have to start making random guesses like suggesting that your r should be different, or that your f3 function should involve x4 instead of x3... because how would we know?
S =
'clc; close all; clear all; S0=100;E0=10;I0=3;R0=0;D0=0; lambda=0;a1=0.3;a2=0.4;r=0.23;K=0.3;phi=0.25;b=0.45;mu=0.12; beta=1;h=0.02;N=5; tfinal=100; f1=@(x1,x2,x3,x4,x5) lambda-a1*x1*x2-a2*x1*x3+r*x4; f2=@(x1,x2,x3,x4,x5) a1*x1*x2+a2*x1*x3-K*x2-phi*x2; f3=@(x1,x2,x3,x4,x5) K*x2-(b+mu)*x3; f4=@(x1,x2,x3,x4,x5) b*x3+phi*x2-r*x4; f5=@(x1,x2,x3,x4,x5) mu*x3; S(N+1)=[0];E(N+1)=[0];I(N+1)=[0];R(N+1)=[0];D(N+1)=[0]; Sp(N+1)=[0];Ep(N+1)=[0];Ip(N+1)=[0];Rp(N+1)=[0];Dp(N+1)=[0]; M1=0;M2=0;M3=0;M4=0;M5=0;L1=0;L2=0;L3=0;L4=0;L5=0; t=0:tfinal/N:tfinal; Sp(1)=S0+h^(beta)*(f1(S0,E0,I0,R0,D0))/gamma(beta)*(beta); Ep(1)=E0+h^(beta)*(f2(S0,E0,I0,R0,D0))/gamma(beta)*(beta); Ip(1)=I0+h^(beta)*(f3(S0,E0,I0,R0,D0))/gamma(beta)*(beta); Rp(1)=R0+h^(beta)*(f4(S0,E0,I0,R0,D0))/gamma(beta)*(beta); Dp(1)=D0+h^(beta)*(f5(S0,E0,I0,R0,D0))/gamma(beta)*(beta); S(1)=S0+h^(beta)*(f1(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f1(S0,E0,I0,R0,D0))/gamma(beta+2); E(1)=E0+h^(beta)*(f2(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f2(S0,E0,I0,R0,D0))/gamma(beta+2); I(1)=I0+h^(beta)*(f3(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f3(S0,E0,I0,R0,D0))/gamma(beta+2); R(1)=R0+h^(beta)*(f4(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f4(S0,E0,I0,R0,D0))/gamma(beta+2); D(1)=D0+h^(beta)*(f5(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f5(S0,E0,I0,R0,D0))/gamma(beta+2); for n=1:N M1=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f1(S0,E0,I0,R0,D0)); M2=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f2(S0,E0,I0,R0,D0)); M3=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f3(S0,E0,I0,R0,D0)); M4=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f4(S0,E0,I0,R0,D0)); M5=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f5(S0,E0,I0,R0,D0)); L1=((n+1)^(beta)-(n)^(beta))*(f1(S0,E0,I0,R0,D0)); L2=((n+1)^(beta)-(n)^(beta))*(f2(S0,E0,I0,R0,D0)); L3=((n+1)^(beta)-(n)^(beta))*(f3(S0,E0,I0,R0,D0)); L4=((n+1)^(beta)-(n)^(beta))*(f4(S0,E0,I0,R0,D0)); L5=((n+1)^(beta)-(n)^(beta))*(f5(S0,E0,I0,R0,D0)); for j=1:n M1=M1+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f1(S(j),E(j),I(j),R(j),D(j))); M2=M2+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f2(S(j),E(j),I(j),R(j),D(j))); M3=M3+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f3(S(j),E(j),I(j),R(j),D(j))); M4=M4+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f4(S(j),E(j),I(j),R(j),D(j))); M5=M5+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f5(S(j),E(j),I(j),R(j),D(j))); L1=L1+((n+1-j)^(beta)-(n-j)^(beta))*(f1(S(j),E(j),I(j),R(j),D(j))); L2=L2+((n+1-j)^(beta)-(n-j)^(beta))*(f2(S(j),E(j),I(j),R(j),D(j))); L3=L3+((n+1-j)^(beta)-(n-j)^(beta))*(f3(S(j),E(j),I(j),R(j),D(j))); L4=L4+((n+1-j)^(beta)-(n-j)^(beta))*(f4(S(j),E(j),I(j),R(j),D(j))); L5=L5+((n+1-j)^(beta)-(n-j)^(beta))*(f5(S(j),E(j),I(j),R(j),D(j))); end Sp(n+1)=(S0+(1/gamma(beta))*L1); Ep(n+1)=(E0+(1/gamma(beta))*L2); Ip(n+1)=(I0+(1/gamma(beta))*L3); Rp(n+1)=(R0+(1/gamma(beta))*L4); Dp(n+1)=(D0+(1/gamma(beta))*L5); S(n+1)=S0+(h)^(beta)*((f1(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M1)/gamma(beta+2); E(n+1)=E0+(h)^(beta)*((f2(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M2)/gamma(beta+2); I(n+1)=I0+(h)^(beta)*((f3(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M3)/gamma(beta+2); R(n+1)=R0+(h)^(beta)*((f4(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M4)/gamma(beta+2); D(n+1)=D0+(h)^(beta)*((f5(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M5)/gamma(beta+2); end Y1(1)=S0; Y2(1)=E0; Y3(1)=I0; Y4(1)=R0; Y5(1)=D0; for k=1:N Y1(k+1)=S(k); Y2(k+1)=E(k); Y3(k+1)=I(k); Y4(k+1)=R(k); Y5(k+1)=D(k); end plot(t,Y1) '
eval(S)

Jan on 3 Dec 2021
Edited: Jan on 3 Dec 2021
I do not know, what you consider as "correct" result. All I see is the running code and there is no chance to predict, what you wnt to do instead.
But at least I can simplify your code, such that a debugging is much easier. Instead of repeating any equation 5 times, I've converted [f1,f2, f3, f4, f5] to one function f and did the same for the 5 inputs:
function Untitled2
lambda=0;a1=0.3;a2=0.4;r=0.23;K=0.3;phi=0.25;b=0.45;mu=0.12;
beta=1;h=0.02;
N = 5;
tfinal = 100;
t = 0:tfinal/N:tfinal;
f = @(x) [ ...
lambda - a1 * x(1) * x(2) - a2 * x(1) * x(3) + r * x(4); ...
a1 * x(1) * x(2) + a2 * x(1) * x(3) - K * x(2) - phi * x(2); ...
K * x(2) - (b + mu) * x(3); ...
b * x(3) + phi * x(2) - r * x(4); ...
mu * x(3)];
S = zeros(5, N+1);
Sp = zeros(5, N+1);
S0 = [100; 10; 3; 0; 0];
Sp(:, 1) = S0 + h^beta * f(S0) / gamma(beta) * beta;
S(:, 1) = S0 + h^beta * f(Sp) / gamma(beta+2) + h^beta * beta *...
f(S0) / gamma(beta+2);
for n=1:N
M = ((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1)) * f(S0);
L = ((n+1)^beta - n^beta) * f(S0);
for j=1:n
M = M + ((n-j+2)^(beta+1) + (n-j)^(beta+1) - 2*(n-j+1)^(beta+1)) * ...
f(S(:, j));
L = L + ((n+1-j)^beta - (n-j)^beta) * f(S(:, j));
end
Sp(:, n+1) = S0 + L / gamma(beta);
S(:, n+1) = S0 + h^(beta) * (f(Sp(:, n+1)) + M) / gamma(beta+2);
end
Y1 = [S0, S(:, 1:N)];
plot(t, Y1)
end
It looks strange, that you calculate S(:, 6), but do not use it in the output.
##### 2 CommentsShowHide 1 older comment
Surath Ghosh on 3 Dec 2021
That is okay. Thank you very much for simplifying. Actually, I am varyfying a research paper. According to that paper values of S will be always less than 100 in the interval (0, 100). Here, I have just plotted S(t) for checking. I have attached the graph. Please check it.

R2016b

### Community Treasure Hunt

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

Start Hunting!

Translated by