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

Question

Surath Ghosh 2021-12-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1602125-please-help-me-to-find-the-errors-in-my-code-for-predictor-corrector-method

评论： Surath Ghosh 2021-12-3

Untitled.m

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2021-12-3

在 MATLAB Online 中打开

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 = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/822230/Untitled.m';

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

eval(S)

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2021-12-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1602125-please-help-me-to-find-the-errors-in-my-code-for-predictor-corrector-method#answer_846695

编辑：Jan 2021-12-3

在 MATLAB Online 中打开

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

Adding spaces around the operators improve the readability.

It looks strange, that you calculate S(:, 6), but do not use it in the output.

2 个评论
显示无隐藏无

Surath Ghosh 2021-12-3

Thank you very much for your comment. Actually Sp, Ip, Ep...... are the predicted values and S, I, E....are the corrected values. I am getting wrong result. I do not understand about the errors in my logic. Can you explain me in details? This is the code of a numerical method that is predictor - Corrector method in mathematics.

Surath Ghosh 2021-12-3

Capture.JPG

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.

请先登录，再进行评论。