当前位置：首页 > news >正文

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

news 来源：原创 2024/9/28 15:35:10

用Adamas内差二步方法，内差三步方法，外差二步方法，外差三步方法这四种方法计算。
在这里插入图片描述

中k为1和2.

在这里插入图片描述
k为2和3

代码

function chap1_adams_methodu0 = 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')endfunction 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:Nf_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_ins2 = 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
endfunction 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:Nf_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_ins2 = 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
endfunction 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
endfunction 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
endfunction f = f1(t, u)
f = -5*u ;
endfunction f = exact1(t)
f = exp(-5*t);
end