% ADAMS VARIABLE STEP-SIZE PREDICTOR-CORRECTOR ALGORITHM 5.5
%
% To approximate the solution of the initial value problem
% y' = f( t, y ), a %
% with local truncation error within a given tolerance:
%
% INPUT: endpoints a, b; initial condition ALPHA; tolerance TOL;
% maximum step size HMAX; minimum step size HMIN.
%
% OUTPUT: I, T(I), W(I), H where at the Ith step W(I) approximates
% y(T(I)) and step size H was used or a message that the
% minimum step size was exceeded.
syms('OK', 'A', 'B', 'ALPHA', 'TOL', 'HMIN', 'HMAX');
syms('FLAG', 'NAME', 'OUP', 'DONE', 'KK', 'NFLAG');
syms('I', 'WP', 'WC', 'SIG', 'K', 'J', 'Q');
syms('W','T','t','y','s','TT','WW','K1','K2','K3','K4');
syms('P1','P2');
TRUE = 1;
FALSE = 0;
% STEP 1 Runge-Kutta Order 4 Method is implemented within the
% following code.
fprintf(1,'This is the Adams Variable Step-size Predictor-');
fprintf(1,'Corrector Method\n');
fprintf(1,'Input the function F(t,y) in terms of t and y\n');
fprintf(1,'For example: y-t^2+1 \n');
s = input(' ','s');
F = inline(s,'t','y');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input left and right endpoints on separate lines.\n');
A = input(' ');
B = input(' ');
if A >= B
fprintf(1,'Left endpoint must be less than right endpoint.\n');
else
OK = TRUE;
end;
end;
OK = FALSE;
fprintf(1,'Input the initial condition.\n');
ALPHA = input(' ');
while OK == FALSE
fprintf(1,'Input tolerance.\n');
TOL = input(' ');
if TOL fprintf(1,'Tolerance must be positive.\n');
else
OK = TRUE;
end;
end;
OK = FALSE;
while OK == FALSE
fprintf(1,'Input minimum and maximum mesh spacing ');
fprintf(1,'on separate lines.\n');
HMIN = input(' ');
HMAX = input(' ');
if HMIN < HMAX & HMIN > 0
OK = TRUE;
else
fprintf(1,'Minimum mesh spacing must be a ');
fprintf(1,'positive real number and less than\n');
fprintf(1,'the maximum mesh spacing.\n');
end;
end;
if OK == TRUE
fprintf(1,'Choice of output method:\n');
fprintf(1,'1. Output to screen\n');
fprintf(1,'2. Output to text file\n');
fprintf(1,'Please enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'For example A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end;
fprintf(OUP, 'ADAMS VARIABLE STEP-SIZE PREDICTOR CORRECTOR METHOD\n\n');
fprintf(OUP, ' t w h sigma\n');
% STEP 2
T = zeros(1,100);
W = zeros(1,100);
T(1) = A;
W(1) = ALPHA;
H = HMAX;
% OK is used in place of FLAG to exit the loop in Step 4.
OK = TRUE;
% DONE is used in place of last to indicate when last value
% is calculated
DONE = FALSE;
% STEP 3
for KK = 1:3
X = T(KK);
Y = W(KK);
K1 = H*F(X,Y);
K2 = H*F(X+0.5*H,Y+0.5*K1);
K3 = H*F(X+0.5*H,Y+0.5*K2);
K4 = H*F(X+H,Y+K3);
WW = Y+(K1+2.0*(K2+K3)+K4)/6.0;
TT = X+H;
T(KK+1) = TT;
W(KK+1) = WW;
end;
% NFLAG indicates the computation from RK4
NFLAG = 1;
I = 5;
% use TT in place of t
TT = T(4) + H;
% STEP 4
while DONE == FALSE
% STEP 5
% predict W(I)
P1 = 55.0*F(T(I-1),W(I-1))-59.0*F(T(I-2),W(I-2));
P2 = 37.0*F(T(I-3),W(I-3))-9.0*F(T(I-4),W(I-4));
WP = W(I-1)+H*(P1+P2)/24.0;
% correct W(I)
P1 = 9.0*F(TT,WP)+19.0*F(T(I-1),W(I-1))-5.0*F(T(I-2),W(I-2));
P2 = F(T(I-3),W(I-3));
WC = W(I-1)+H*(P1+P2)/24.0;
SIG = 19.0*abs(WC-WP)/(270.0*H);
% STEP 6
if SIG % STEP 7
% result accepted
W(I) = WC;
T(I) = TT;
% STEP 8
if NFLAG == 1
K = I-3;
KK = I-1;
% Previous results are also accepted.
for J = K:KK
fprintf(OUP, '%12.8f %11.8f %11.8f %11.8f\n', T(J), W(J), H, SIG);
end;
fprintf(OUP, '%12.8f %11.8f %11.8f %11.8f\n', T(I), W(I), H, SIG);
else
% Previous results were already accepted.
fprintf(OUP, '%12.8f %11.8f %11.8f %11.8f\n', T(I), W(I), H, SIG);
end;
% STEP 9
if OK == FALSE
% Next step is 20.
DONE = TRUE;
else
% STEP 10
I = I+1;
NFLAG = 0;
% STEP 11
if SIG B
% Increase H if more accuracy than required has been obtained,
% or decrease H to include b as a mesh point.
% STEP 12
% to avoid underflow
if SIG Q = 4.0;
else
Q = (0.5*TOL/SIG)^(1/4);
end;
% STEP 13
if Q > 4.0
H = 4.0*H;
else
H = Q * H;
end;
% STEP 14
if H > HMAX
H = HMAX;
end;
% STEP 15
if T(I-1)+4.0*H > B
H = 0.25*(B-T(I-1));
if H < TOL
DONE = TRUE;
end;
OK = FALSE;
end;
% STEP 16
for KK = I-1:I+2
X = T(KK);
Y = W(KK);
K1 = H*F(X,Y);
K2 = H*F(X+0.5*H,Y+0.5*K1);
K3 = H*F(X+0.5*H,Y+0.5*K2);
K4 = H*F(X+H,Y+K3);
WW = Y+(K1+2.0*(K2+K3)+K4)/6.0;
TT = X+H;
T(KK+1) = TT;
W(KK+1) = WW;
end;
NFLAG = 1;
I = I+3;
end;
end;
else
% FALSE branch for Step 6 - result rejected.
% STEP 17
Q = (0.5*TOL/SIG)^(1/4);
% STEP 18
if Q < 0.1
H = 0.1 * H;
else
H = Q * H;
end;
% STEP 19
if H < HMIN
fprintf(OUP, 'HMIN exceeded\n');
DONE = TRUE;
else
if T(I-1)+4.0*H > B
H = 0.25*(B-T(I-1));
end;
if NFLAG == 1
% Previous results also rejected.
I = I-3;
end;
for KK = I-1:I+2
X = T(KK);
Y = W(KK);
K1 = H*F(X,Y);
K2 = H*F(X+0.5*H,Y+0.5*K1);
K3 = H*F(X+0.5*H,Y+0.5*K2);
K4 = H*F(X+H,Y+K3);
WW = Y+(K1+2.0*(K2+K3)+K4)/6.0;
TT = X+H;
T(KK+1) = TT;
W(KK+1) = WW;
end;
I = I+3;
NFLAG = 1;
end;
end;
% STEP 20
TT = T(I-1) + H;
end;
% STEP 21
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully \n',NAME);
end;
end;
