Adams方法：

matlab代码：

function chap1_adams_method
% test the Adams interpolation method and the Adams extrapolation method for ODE IVP
% /   u' = f(t,u),
% \   u(0) = u0.

% foundate = '2015-3-8';
% chgedate = '2020-04-22';
% by Zhang, Xuping

u0 = 1;
T = 2;
h = 0.1;
N = T/h;
t = 0:h:T;
solu = exact1(t);

f = @f1;
u_inter_2s = adams_inter_2steps(f,u0,t,h,N);
u_extra_2s = adams_extra_2steps(f,u0,t,h,N);
figure(1)
plot(t,u_inter_2s,'*',t,u_extra_2s,'o',t,solu,'r')
legend('Adams-inter-2s','Adams-extra-2s','Exact-soln')

u_inter_3s = adams_inter_3steps(f,u0,t,h,N);
u_extra_3s = adams_extra_3steps(f,u0,t,h,N);
figure(2)
plot(t,u_inter_3s,'*',t,u_extra_3s,'o',t,solu,'r')
legend('Adams-inter-3s','Adams-extra-3s','Exact-soln')

end
% ----------------------------- subroutines -------------------------------
function u = adams_inter_2steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
u(2) = exact1(1*h);
eps_in = 1e-6;
K_in = 6;
for n = 2:N
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    s1 = u(n);
    du = 1;
    k = 1;
    while abs(du)>eps_in & k<K_in
        s2 = u(n) + h*( 5*f(t(n+1),s1) + 8*f_n - f_nm1 )/12;
        du = s2 - s1;
        s1 = s2;
        k = k + 1;
    end
    u(n+1) = s2;
end
end

function u = adams_inter_3steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
% u(3) = u(2) + h*f(t(2),u(2));
u(2) = exact1(1*h);
u(3) = exact1(2*h);
eps_in = 1e-6;
K_in = 6;
for n = 3:N
    f_nm2 = f(t(n-2),u(n-2));
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    s1 = u(n);
    du = 1;
    k = 1;
    while abs(du)>eps_in & k<K_in
        s2 = u(n) + h*( 9*f(t(n+1),s1) + 19*f_n - 5*f_nm1 + f_nm2 )/24;
        du = s2 - s1;
        s1 = s2;
        k = k + 1;
    end
    u(n+1) = s2;
end
end

function u = adams_extra_2steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
u(2) = exact1(h);
for n = 2:N
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    u(n+1) = u(n) + h*( 3*f_n - f_nm1 )/2;
end
end

function u = adams_extra_3steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
% u(3) = u(2) + h*f(t(2),u(2));
u(2) = exact1(1*h);
u(3) = exact1(2*h);
eps_in = 1e-6;
for n = 3:N
    f_nm2 = f(t(n-2),u(n-2));
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    u(n+1) = u(n) + h*( 23*f_n - 16*f_nm1 + 5*f_nm2 )/12;
end
end

function f = f1(t,u)
f = -5*u;
end

function f = exact1(t)
f = exp(-5*t);
end

function chap1_adams_method_stability
% test the stability of adams method
% /   u' = f(t,u),
% \   u(0) = u0.

u0 = 1;
T = 5;
h = 0.01;
N = T/h;
t = 0:h:T;
solu = exact1(t);

f = @f1;
u_inter_3s = adams_inter_3steps(f,u0,t,h,N);
u_extra_4s = adams_extra_4steps(f,u0,t,h,N);
figure(1)
plot(t,u_inter_3s,'*',t,u_extra_4s,'o',t,solu,'r')
legend('Adams-inter-3s','Adams-extra-4s','Exact-soln')

end
% ----------------------------- subroutines -------------------------------
function u = adams_inter_3steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
% u(3) = u(2) + h*f(t(2),u(2));
u(2) = exact1(1*h);
u(3) = exact1(2*h);
eps_in = 1e-6;
K_in = 6;
for n = 3:N
    f_nm2 = f(t(n-2),u(n-2));
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    s1 = u(n);
    du = 1;
    k = 1;
    while abs(du)>eps_in & k<K_in
        s2 = u(n) + h*( 9*f(t(n+1),s1) + 19*f_n - 5*f_nm1 + f_nm2 )/24;
        du = s2 - s1;
        s1 = s2;
        k = k + 1;
    end
    u(n+1) = s2;
end
end

function u = adams_extra_2steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
u(2) = exact1(h);
for n = 2:N
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    u(n+1) = u(n) + h*( 3*f_n - f_nm1 )/2;
end
end

function u = adams_extra_3steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
% u(3) = u(2) + h*f(t(2),u(2));
u(2) = exact1(1*h);
u(3) = exact1(2*h);
eps_in = 1e-6;
for n = 3:N
    f_nm2 = f(t(n-2),u(n-2));
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    u(n+1) = u(n) + h*( 23*f_n - 16*f_nm1 + 5*f_nm2 )/12;
end
end

function u = adams_extra_4steps(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
% u(2) = u(1) + h*f(t(1),u(1));
% u(3) = u(2) + h*f(t(2),u(2));
% u(4) = u(3) + h*f(t(3),u(3));
u(2) = exact1(h);
u(3) = exact1(2*h);
u(4) = exact1(3*h);
for n = 4:N
    f_nm3 = f(t(n-3),u(n-3));
    f_nm2 = f(t(n-2),u(n-2));
    f_nm1 = f(t(n-1),u(n-1));
    f_n = f(t(n),u(n));
    u(n+1) = u(n) + h*( 55*f_n - 59*f_nm1 + 37*f_nm2 - 9*f_nm3 )/24;
end
end

function f = f1(t,u)
f = -5*u;
end

function f = exact1(t)
f = exp(-5*t);
end

注：课件内容为大连理工大学张旭平老师所写