matlab线性多部法求常微分方程数值解

用Adamas内差二步方法，内差三步方法，外差二步方法，外差三步方法这四种方法计算。


中k为1和2.

k为2和3

代码

function chap1_adams_method

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

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

结果：