How i implement Adams Predictor-Corrector Method from general code ?

55 次查看(过去 30 天)
Below is the Adams predictor-corrector formula and general code. How can I adapt this code to the above question? Can you please help?
%------------------------------------------------------
% 2-step Predictor-Corrector
% [T,Y]=dd2(f,definition,y,h); definition=[t1,tfinal]
%------------------------------------------------------
function [T,Y]=dd2(f,definition,Y1,h)
t1=definition(1);tfinal=definition(2);T=t1;Y=Y1;
t2=t1+h;
definition=[t1,t2];
[T,Y]=rk2(f,definition,Y1,h) ;
Y2=Y(2);
while t2 <tfinal
t3=t2+h;
P=Y2+h*(3/2*f(t2,Y2)-1/2*f(t1,Y1));
Y3=Y2+h/12*(5* f(t3,P)+8*f(t2,Y2)-f(t1,Y1));
Y1=Y2; Y2=Y3;t1=t2;t2=t3;
T=[T;t3];Y=[Y;Y3];
end
%
%----------------------------------------------
  2 个评论
Torsten
Torsten 2022-5-27
You know the correct result of your differential equation.
If you plot Y against T in the calling program and compare the plot with the analytical solution, both should be approximately the same.
If yes, your code is (most probably) correct, if not, it's not.

请先登录,再进行评论。

采纳的回答

Lateef Adewale Kareem
编辑:Lateef Adewale Kareem 2022-5-30
clc; clear all;
h = 0.01;
mu = 20;
f_m = @(t,y) mu*(y-cos(t))-sin(t);
exact = @(t) exp(mu*t)+cos(t);
[t,y_m] = dd2(f_m,[0, 1],exact(0), exact(h), h);
plot(t, exact(t)); hold
Current plot held
plot(t,y_m);
%plot(t,y,'-o');
legend('Exact Solution','Adams predictor-corrector formula')
xlabel('t')
ylabel('y')
title('When h = 0.01 and µ=20')
%------------------------------------------------------
% 2-step Predictor-Corrector
% [T,Y]=dd2(f,definition,y,h); definition=[t1,tfinal]
%------------------------------------------------------
function [T,Y] = dd2(f, definition, Y1, Y2, h)
t1 = definition(1); tfinal = definition(2); t = t1:h:tfinal;
T = t(1:2)'; Y = [Y1;Y2];
for i = 2:numel(t)-1
P = Y(i) + h/2*(3*f(t(i),Y(i))-f(t(i-1),Y(i-1)));
Y(i+1) = Y(i) + h/12*(5*f(t(i+1), P) + 8*f(t(i),Y(i)) - f(t(i-1),Y(i-1)));
T=[T;t(i+1)];
end
end
%
  4 个评论
Lateef Adewale Kareem
编辑:Lateef Adewale Kareem 2022-5-30
yeah. he should have sent it in. I have modified the solution to use the exact solution at h

请先登录,再进行评论。

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by