[微]显式Euler方法，梯形方法，预估-校正Euler方法参考程序及解释

程序中的注释充分详细。大部分同学没有解决好 梯形方法中求隐式方程的部分，请参考程序中的while循环。改程序给出了5种不同步长下的求解区间右端点处的误差值，分别存在了errF,errT,errPC中，程序运行之后会给出errF,errT,errPC及5幅对比图像。通过errF中误差随步长减半而减半可知显式Euler方法整体截断误差阶为1；errT和errPC中的数值表明，步长减半时误差逐减为原1/4，可知梯形方法和预估校正Euler方法的整体截断误差阶为2。

代码:

% 本程序实现求解下述问题的 显式Euler方法，梯形Euler方法，预估校正Euler方法
% 一阶常微分方程初值问题： y' = - y^2,  y(0)=1, [0,1], 精确解y(x)=1/(1+x)

clear;clc;close all;
format long
% errF,errT,errPC分别记录三种方法在5种步长下 求解区间右端点处的误差
errF = zeros(5,1);
errT = zeros(5,1);
errPC = zeros(5,1);
n = 8;
for M = 1:5
    % 用下式来给出不同的步长h=1/6,1/32,1/64,1/128,1/256
    n = n*2;
    h = 1/n;
    x = 0:h:1;
       
    % 精确解
    y_exact = zeros(n+1,1);
    for i = 1:n+1
        y_exact(i) = 1/(1+x(i));
    end

    % 显式Euler方法
    yF = zeros(n+1,1);
    yF(1) = 1;
    for i = 1:n
        yF(i+1) = yF(i) - h*yF(i)^2;
    end

    % 梯形方法
    yT = zeros(n+1,1);
    yT(1) = 1;
    for i = 1:n
        % 以这种方式给初值
        yT(i+1) = yT(i) - h*yT(i)^2;
        % 简单迭代求解隐式方程
        err = 1;
        k = 0;
        while err> 1.0e-9
            temp = yT(i) - h*(yT(i+1)^2 + yT(i)^2)/2;
            err = abs(temp - yT(i+1));
            yT(i+1) = temp;
            % 用k记录迭代次数
            k = k + 1;
        end
    end

    % 预估-校正Euler方法
    yPC = zeros(n+1,1);
    yPC(1) = 1;
    for i = 1:n
        yPC(i+1) = yPC(i) - h*yPC(i)^2;
        yPC(i+1) = yPC(i) - h*(yPC(i+1)^2 + yPC(i)^2)/2;
    end

    figure
    plot(x,y_exact,'r',x,yF,'-.',x,yT,'+',x,yPC,'*')
    legend('精确','显式Euler','梯形','预校Euler')
   
   
    errF(M) = abs(yF(n+1) - y_exact(n+1));
    errT(M) = abs(yT(n+1) - y_exact(n+1));
    errPC(M) = abs(yPC(n+1) - y_exact(n+1));
end
errF
errT
errPC