【微分方程数值解】四阶Adams公式的显式隐式的实现与比较

·实验原理

 ·实验步骤

 
 分析与讨论

1.收敛性：当步长h取得充分小时，数值解能否足够精确地逼近初值问题的真解

2.局部截断误差：通过前面求解的节点结果推出后续节点的值，局部截断误差为则称该算法具有p阶精度

3.稳定性：最初产生的误差在以后各步的计算中是否会无限制扩大

思考与回答

1.在h变小的过程中，四阶显式Adams算法和隐式Adams算法的精度逐渐升高，接近于精确解。

2.四阶显式隐式Adams局部截断误差：

3.隐式方法的精度及稳定性比显式方法好得多

代码实现

function Adam
T = 1;
h = 0.1;
t = 0:h:T;
N = length(t)-1;
solu = exp(-5.0*t);
u0 = 1;
f = @f1;
u_out_Adams = Adams1(f,u0,t,h,N);
u_in_Adams = Adams2(f,u0,t,h,N);
 
figure (1)
plot(t,u_out_Adams,'*',t,solu,t,u_in_Adams,'d');
legend('显示Adams','精确解','隐式Adams');
end
 
function u = euler(f,u0,t,h,N)
u    = zeros(N+1,1);
u(1) = u0;
for n = 1:N
    fn     = f(t(n),u(n));
    u(n+1) = u(n)+h*fn;
end
end
function u = Adams1(f,u0,t,h,N)
u=zeros(N+1,1);
u(1)=u0;
u = euler(f,u0,t,h,4);
for n=4:N
    f1 = f(t(n),u(n));
    f2 = f(t(n-1),u(n-1));
    f3 = f(t(n-2),u(n-2));
    f4 = f(t(n-3),u(n-3));
    u(n+1) = u(n)+(h/24)*(55*f1-59*f2+37*f3-9*f4);
end
end
function u = Adams2(f, u0, t, h ,N)
u    = zeros(N+1,1);
u(1) = u0;
u = euler(f,u0,t,h,4);
eps_in = 1e-6;
  K_in = 6;
for n = 3:N
    s1 = u(n);
    du = 1;
    k  = 1; 
    f2 = f(t(n),u(n));
    f3 = f(t(n-1),u(n-1));
    f4 = f(t(n-2),u(n-2));
    while abs(du)>eps_in && k<K_in
        f1 = f(t(n+1),s1);
        s2 = u(n)+(h/24)*(9*f1+19*f2-5*f3+f4);
        du = s2-s1;
        s1 = s2;
         k = k+1; end
    u(n+1) = s2;
end
end
function f =f1(t,u)
f = -5*u;
end

 注：之后将h分别改为h=0.05、h=0.01，得出步骤中的三个图像进行收敛性的比较