Please help me to find the errors in my code for predictor corrector method.

3 views (last 30 days)
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?
  1 Comment
Walter Roberson
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?
filename = '';
S = urlread(filename)
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) '

Sign in to comment.

Answers (1)

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
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));
Sp(:, n+1) = S0 + L / gamma(beta);
S(:, n+1) = S0 + h^(beta) * (f(Sp(:, n+1)) + M) / gamma(beta+2);
Y1 = [S0, S(:, 1:N)];
plot(t, Y1)
Adding spaces around the operators improve the readability.
It looks strange, that you calculate S(:, 6), but do not use it in the output.
Surath Ghosh
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.

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!