Francesco Tudisco

Associate Professor (Reader) in Machine Learning

School of Mathematics, The University of Edinburgh
The Maxwell Institute for Mathematical Sciences
School of Mathematics, Gran Sasso Science Institute JCMB, King’s Buildings, Edinburgh EH93FD UK
email: f dot tudisco at ed.ac.uk

Matlab example 2: Explicit vs Implicit 2 steps Adams methods

I recommend you try to write a script that implements Adams' methods with an arbitrary number of steps. Then use this to make a numerical evaluation of what you observed with homework no. 1.

clear all
close all

f = @(x,y) y^2;
real_y = @(x) 1/(1-x);

I = [0 0.99];
numsteps = 20;
h = (I(2)-I(1))/numsteps;

x(1) = 0;
y_e(1) = 1;
y_i(1) = 1;
y(1) = 1;

figure
plot(x,y_e,'-*b');
hold on
plot(x,y_i,'-*r');
plot(x,y,'--k');
xlim(I)
ylim([0, 20])
axis square

F_e = [f(x(1),y_e(1))];
F_i = [f(x(1),y_i(1))];

for i = 1 : numsteps-1
    % Compute the new grid point
    x(i+1) = x(i) + h;
    % Compute the new approximation with EAM
    y_e(i+1) = EAM2(y_e(i),F_e,h);
    F_e = [f(x(i+1),y_e(i+1)) F_e];
    % Compute the new approximation with IAM
    y_i(i+1) = PEC_IAM2(x(i+1),y_i(i),F_i,h,f);
    F_i = [f(x(i+1),y_i(i+1)) F_i];
    % Compute the value of the true solution
    y(i+1) = real_y(x(i+1));

    plot(x,y_e,'-*b');
    plot(x,y_i,'-*r');
    plot(x,y,'--k');
    xlim(I)
    ylim([0, 20])
    axis square
    pause(.1)
end
legend('EAM(2)', 'IAM(2) with PEC', 'True Solution', 'location','northwest');

hold off






function y = EAM2(yn, F, h)

    % INPUT:   yn : solution approx on x_n
    %          F  : vector of previous values F = [f_n, f_{n-1}]
    %          h  : step size
    % OUTPUT:  y  : new value y_{n+1}

    if length(F) >= 2
        y = yn + h * .5 * (3*F(1) - F(2));
    elseif length(F) == 1
        y = yn + h * F(1); %Explicit Euler method
    end    

end

function y = PEC_IAM2(xn1, yn, F, h, f)

    % INPUT:   xn1 : next grid point x_{n+1}
    %          yn  : solution approx on x_n
    %          F   : vector of previous values F = [f_n, f_{n-1}]
    %          h   : step size
    %          f   : function handler for f(x,y)
    % OUTPUT:  y   : new value y_{n+1}

    % Predictor
    yn1 = EAM2(yn, F, h);

    % Evaluation
    F = [f(xn1,yn1) F];

    % Corrector (IAM2)
    y = yn + h * .5 * (F(1) + F(2));

end