Variable step size Implicit Adams method
Using the recursions for the coefficients Φj(n), Ψj(n) and gj(n), which we recall below
Φ0(n)=Ψ0(n)=fn
Φj+1(n)=Φj(n)−Ψj(n−1)
Ψj+1(n)=βj+1(n)Φj+1(n)
βj(n)=βj−1(n)(xn+1−xn−j+1xn−xn−j), β0(n)=1
Q(0,q)=1/q, Q(1,q)=1q(q+1)
Q(j,q)=Q(j−1,q)−Q(j−1,q+1)
gj(n)=Q(j,1)
Implement a 1-step and 2-step IAM with adaptive step size to solve the ODE y′=y2, y(0)=1, using the scheme
yn+1=pn+1+hngk(n)Φk(n+1), pn+1=yn+hn∑k−1j=0gj(n)Ψj(n)
Matlab DEMO
clear all
% F = history of function evaluations
% X = history of points
% Y = history of solution approximations
f = @(x,y) y^2;
sol = @(x) 1./(1-x); %true solution
X = [0.01 0.02 0.03];
Y = [sol(X(1)) sol(X(2)) sol(X(3))];
F = [f(X(1),Y(1)) f(X(2),Y(2)) f(X(3),Y(3))];
hnew = 0.13;
h = hnew;
tol = 1e-5;
% Interval [0,1]
while X(end) < 0.99
[y,hnew] = corrector(hnew,F,X,Y,f,tol);
X = [X X(end) + hnew];
Y = [Y y];
F = [F f(X(end),Y(end))];
h = [h hnew];
end
Y = Y(1:end-1);
X = X(1:end-1);
close all
subplot(121)
plot(X,Y,'-or');
hold on
plot(X,sol(X),'--k', 'linewidth',3)
title('Approx solution')
subplot(122)
semilogx(h,'-')
title('Sequence of steps h_n')
function [y,hnew] = corrector(hnew,F,X,Y,f,tol)
for it = 1:100
[p,PSI,G] = predictor(hnew,F,X,Y);
Xn1 = X(end)+hnew;
G(3) = 1/2 - (1/6)*(hnew/(Xn1-X(end-1)));
PHI3 = f(Xn1,p) - F(end) - PSI(2);
y = p + hnew * G(3) * PHI3;
hn = hnew;
%%% We compute an approximation LE for the local error at the next point
%%% and compute a new step length as a consequence
% compute PHI_3(n+1)
b1n = hn / (X(end)-X(end-1));
c1 = (Xn1 - X(end-1))/(X(end)-X(end-2));
c2 = (X(end)-X(end-1))/(X(end-1)-X(end-2));
PHI4 = PHI3 - b1n * c1 * (F(end)-F(end-1)-c2*(F(end-1)-F(end-2)));
% compute g_3(n)
G(4) = G(3) - (hn/(Xn1-X(end-2))) * (1/6 - 1/12 * (hn / (Xn1 - X(end-1))));
% compute LE
LE = hn * abs(PHI4 * (G(4)-G(3)) );
% define a new hnew
hnew = hn * (tol/LE)^(1/4);
%%% If the approximate LE is smaller than the tollerance we accept
if LE < tol
break;
end
end
end
function [p,PSI,G] = predictor(hnew,F,X,Y)
G(1) = 1;
G(2) = 1/2;
PSI(1) = F(end);
Xn1 = X(end) + hnew;
PSI(2) = (F(end)-F(end-1))*((Xn1-X(end))/(X(end)-X(end-1)));
p = Y(end) + hnew * (G(1)*PSI(1) + G(2)*PSI(2));
end