数值分析常见基本算法及MATLAB代码总结

刚考完研究生的数值分析课，打算整理一下平时的上机代码，做一个汇总。

数值分析大致八部分内容——

解线性方程组的直接法：Gauss消去法与矩阵三角分解法（Doolittle分解法相比Crout分解法更常用）及其选择列主元的改进方法、Doolittle分解法的延伸（实对称正定矩阵利用Cholesky分解得到的平方根法、三对角矩阵作为线性方程组系数矩阵的追赶法）

解线性方程组的迭代法：Jacobi迭代法、Gauss-Seidel迭代法（利用前者每次迭代已得到的最新分量加速）、逐次超松弛（SOR，Successive Over-Relaxation）方法（取适当松弛因子ω可以比Gauss-Seidel方法收敛得更快）以及共轭梯度法（教材未讲，前提条件：系数矩阵A为对称正定矩阵）

函数拟合的插值法：拉格朗日（Lagrange）插值法与牛顿（Newton）插值法（根据插值多项式的存在唯一性定理，对同一插值问题两者展开的结果相同。牛顿插值法的优势在于：每次增加或删减一个节点，只需插值多项式增加或减少一项，相当于对不同插值节点数的情形给出了迭代公式）、分段插值法（克服插值节点过多时逼近效果反而很差的Runge现象，常见的是分段线性插值，也包括分段抛物线插值法）、Hermite插值法（插值条件包括所有已知节点的函数值和导数值，即两种条件的等式数量相等）、样条插值法（最基本的是三次样条插值法，也属于分段插值法的一种，整体具有二阶连续微商）。需要指出的是，除了分段插值法和样条插值法，其他几种方法经常是交叉混合使用的，因为实际插值条件中节点函数值和节点导数值条件的数目往往不同，可以灵活选择部分条件先构造Lagrange或Hermite插值多项式，然后再添加一项，使得满足附加插值条件的同时，已有插值条件不受影响。两项相加并解出待定系数就是完整的一般形式的插值多项式。此时插值余项也需自行推导，但一般是反复运用罗尔定理得出的带f(ξ)的高阶导数的表达式。

函数近似替代（复杂化简单、高次化低次）的数值逼近方法：数值逼近中引入了函数范数和函数内积的概念。前者用来度量逼近函数与原函数在一个区间内的整体误差，后者广泛用于各种数值逼近方法的计算过程中。函数的∞-范数对应最佳一致逼近；函数的2-范数（Euclid-范数）对应最佳平方逼近。这两种是最基本的数值逼近方法。（上课只选讲了最佳平方逼近，按老师原话，对工科学生来说最佳一致逼近理论性太强，且实际中用得不是很多）。将函数逼近中的线性无关函数族（类比线性代数中n维向量空间的向量线性表示选择的一组基）进行正交化（类似于线性代数中的施密特正交化），取适当的积分区间和权函数，可以衍生出最佳一致逼近中的切比雪夫（Chebyshev）多项式和最佳平方逼近中的勒让德（Legendre）多项式。以引入勒让德多项式为例，此时解法方程组，系数矩阵G为对角阵，可以大大降低函数逼近计算中解线性方程组的工作量。使用的教材中，将统计学中数据拟合的最小二乘法也归到数值逼近一章，个人理解是因为它的理论证明过程和计算套路都和最佳平方逼近及其相似，可以类比学习。

数值积分算法与数值微分（未学）：数值积分一章从机械求积公式出发，引入插值法并采用等距插值节点得到Newton-Cote公式。积分区间等分数n取不同整数有一系列不同数值求积公式，常见的为梯形公式（n=1)、Simpson公式（n=2）及Cotes公式（n=4）。当n≥8时，Cotes系数出现负数，数值不稳定，又有与Newton-Cotes公式对应的一系列复化求积公式（为保证数值积分精度，需积分步长相对固定，维持在一个较小的值。可以等分积分区间，在每个子区间中用n较小(＜8)的Newton-Cotes公式，这是复化求积公式的基本思想。）除此之外数值积分还有变步长的Romberg方法，适用于无法事先选择恰当的步长，在计算过程需要改变步长的情形。

解非线性方程及方程组的数值方法：对于非线性方程，最简单的方法是二分法，这个在高中就涉及过。常用的迭代法包括不动点迭代法、牛顿迭代法（Newton-Raphson迭代法或N-R法）和割线法。对于不动点迭代法，需将f(x)=0适当变形，取合适的φ（x）满足不动点迭代法的收敛条件（①"压缩镜像“；②迭代初值选取区间上导数值的绝对值恒不大于一个小于1的正数）；N-R法可以看做不动点迭代法φ（x）取一特殊函数的特例，也具有鲜明的几何意义，又称切线法；割线法是考虑到部分函数计算导数困难，用差商近似导数所做的改进。Newton法可以推广到非线性方程组的情形，在此基础上有一类拟Newton法，是针对前者的Jacobian矩阵中精确求解数目为方阵阶数n的平方的偏导数值的困难而做的改进。

矩阵特征值的数值解法：乘幂法与反幂法、HouseHolder方法、QR分解法（上课只讲了乘幂法，其他的以后用到了再学习。）

常微分方程的数值解法：欧拉（Euler）方法(前进与后退的Euler公式，综合两者的梯形公式，以及梯形公式中u_k+1项再代入前进的Euler公式，消除隐式参数，得到改进的欧拉公式)、龙格-库塔（Runge-Kutta）法、线性多步法。（P.S. 龙格-库塔法虽然是选学内容，上课没讲，但是对于电气、计算机、土木等工科研究生来说可能是一个重要的数学工具，给他个triple star 吧，今后用到了会细究。印象中本科的电器学教材中也有一节提到了这个）

数值分析课程中还有一些进阶模块，例如积分方程或积分-微分方程的数值解法，各种最优化方法等，也许超出了研究生教学大纲要求掌握的基本内容，今后如有用到再做补充。

下面逐一上代码，有些数值方法可能理论性较强，教材没有给出示例MATLAB程序或还没有MATLAB运行验证过，日后再更新。

1、选取列主元的高斯消去法

function x= pivot_Gauss(A,b)
%功能：用Gauss列主元消去法解n阶线性方程组Ax=b
n=length(b);
for k=1:n-1 %步骤2：消元
    %步骤2.1：选列主元
    [max_value,max_index]=max(abs(A(k:n,k)));
    rk = k+max_index-1; %加上该列中前n-1个元素，还原主元在矩阵A中的行标
    if max_value == 0 %步骤2.2：若A(rk,k)=0,则输出信息“A为奇异阵”，算法终止；否则执行步骤2.3
        warning('系数矩阵奇异！');
        return;
    end
    if rk~=k %步骤2.3：若rk!=k，则交换增广矩阵的第k行和第rk行
        t=A(k,:);A(k,:)=A(rk,:);A(rk,:)=t;
        %交换矩阵的第k行与第rk行（从第k列开始即可，因为前k-1列下三角元素均已消零）
        t=b(k);b(k)=b(rk);b(rk)=t;
    end
    for i=k+1:n %步骤2.4：对i=k+1,...,n做以下两步：计算乘子l_ik;计算a_ij和b_i
        L(i,k)=A(i,k)/A(k,k);
        A(i,k+1:n)=A(i,k+1:n)-L(i,k)*A(k,k+1:n);
        b(i)=b(i)-L(i,k)*b(k);
    end
end
%步骤3：若a_nn=0，则输出信息“A为奇异阵！”，算法终止
if A(n,n)==0
    warning('系数矩阵奇异！');
    return;
end
%步骤4：回代求解
for k=n:-1:1
    if k==n
        x(n)=b(n)/A(n,n);
    else
        x(k)=(b(k)-sum(A(k,k+1:n).*x(k+1:n)))/A(k,k);
    end
end
x

2、选取列主元的Doolittle分解法

% Doolittle.m
function x=Doolittle2_3(A,b)
%功能：用按列选主元的Doolittle方法求解n阶线性方程组Ax=b
%假设A的各阶顺序主子式均不为零
n=length(A);L=eye(n);U=zeros(n);
%length(A)返回矩阵A最大维数；
%I = eye(n) 返回一个主对角线元素为 1 且其他位置元素为 0 的 n×n 单位矩阵。

%步骤1：输入系数矩阵A和及右端向量b
%步骤2：对k=1,2,...,n-1执行步骤2.1~2.4
for k=1:n-1 %省略了中间步长为1
    for i=k:n
        %步骤2.1 计算下三角矩阵L第k列元素的一部分
        s(i)=A(i,k)-L(i,1:k-1)*U(1:k-1,k);
    end
    %步骤2.2 选列主元
    [s_q,q]=max(abs(s(k:n)));
    %[M,I] = max(___) 使用先前语法中的任何输入参数，查找 A 的最大值的索引，并在输出矢量 I 中返回这些索引。
    %如果最大值出现多次，则 max 返回对应于第一次出现位置的索引。
    q=q+k-1;%还原s(i)的索引下标
    %步骤2.3 若q!=k
    if q>k %绝对值最大的列主元不在第q行即当前第一行
    %步骤2.3.1：交换A的第k行与第q行
        t_A=A(k,:);A(k,:)=A(q,:);A(q,:)=t_A;
    %步骤2.3.2：交换下三角矩阵L第k列元素s_k与s_q（l_kk与l_qk）
        t_s=s(k);s(k)=s(q);s(q)=t_s;
    %步骤2.3.3：将下三角矩阵L的第k行前(k-1)个元素与第q行的前(k-1)个元素交换
        t_L=L(k,1:k-1);L(k,1:k-1)=L(q,1:k-1);L(q,1:k-1)=t_L;
    %步骤2.3.4：交换b_k与b_q
        t_b=b(k);b(k)=b(q);b(q)=t_b;
    end
    %步骤2.4：计算矩阵U的第k行元素及矩阵L的第k列元素
    U(k,k)=s(k);
    for i=k+1:n
        U(k,i)=A(k,i)-L(k,1:k-1)*U(1:k-1,i);
        L(i,k)=s(i)/U(k,k)
    end
end
%步骤3：求解上三角矩阵第n行的唯一待定元素u_nn
U(n,n)=A(n,n)-L(n,1:n-1)*U(1:n-1,n);
%步骤4：解单位下三角方程组Ly=b
y=zeros(n,1);
for k=1:n
    if k==1
        y(1)=b(1);
    else 
        y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
    end
end
%步骤5：解上三角方程组Ux=y;
x=zeros(n,1);
for k=n:-1:1
    if k==n
        x(n)=y(n)/U(n,n);
    else
        x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
    end
end
%步骤6：输出解x
x=x';%把列向量转置成行向量，尽量不占用命令行的显示空间。

3、解线性方程组系数矩阵为三对角矩阵的追赶法

%ForwardBackward.m 追赶法解满足行对角占优条件的三对角矩阵的MATLAB程序
%行对角占优条件：(1)比较绝对值，第一行的主对角元大于第一行的另一个非零元，且绝对值均不为零；
%(2)比较绝对值，最后一行的主对角元大于最后一行的另一个非零元，且绝对值均不为零；
%(3)第一行和最后一行之外的其余各行，主对角元的绝对值大于另外两个非零元的绝对值之和且均不为0。

%说明：矩阵A的主对角线元素为d_1~d_n;次下对角线元素为a_2~a_n;次上对角线元素为c_1~c_n-1。

function x=ForwardBackward2_3(a,b,c,d)
n=length(d);
%步骤1：输入4个一维数组a,c,d(存放A的三条对角线上的元素),b(Ax=b，列向量)，由算例程序给出
%步骤2：若d_1==0，则输出信息“方法失败！”，算法终止。
if d(1)==0
    warning('方法失败！')
    return
end
%步骤3：计算p1=d1,q1=c1/d1
p(1)=d(1);q(1)=c(1)/d(1);
%步骤4：对i=2,3,...,n-1，执行步骤4.1~4.3
for i=2:n-1
    %步骤4.1：计算p_i
    p(i)=d(i)-a(i)*q(i-1);
    %步骤4.2：若p_i=0,则输出信息“方法失败！"，算法终止(可以证明系数矩阵A满足对角占优条件，则p_1~p_n均不为0)
    if p(i)==0
        warning('方法失败！')
    return
    end
    %步骤4.3：计算q_i
    q(i)=c(i)/p(i);
end
%步骤5：计算p_n=d_n-a_n*q_n-1。
p(n)=d(n)-a(n)*q(n-1);
%步骤6：若p_n=0，则输出信息“方法失败！”，算法终止。
if p(n)==0
    warning('方法失败！')
    return
end
%步骤7：解下三角方程组Ly=b
y(1)=b(1)/p(1);
for i=2:n
    y(i)=(b(i)-a(i)*y(i-1))/p(i);
end
%步骤8：解上三角方程组Ux=y
x(n)=y(n);
for i=n-1:-1:1
    x(i)=y(i)-q(i)*x(i+1);
end
x=x'

4、Jacobi迭代法与 Gauss-Seidel迭代法

function x = Jacobi3_2(A,b,x0,eps,N)
%功能：用Jacobi迭代法解n阶线性方程组Ax=b
n=length(b);x=ones(n,1);k=0;
%步骤1：输入系数矩阵A，右端向量b,以及初始向量x0，精读eps,以及最大迭代次数N；令k=1

%步骤2：当k≤N时，执行步骤2.1~2.3。
while k<=N
    %步骤2.1：算出x1~xn的第k次迭代递推式
    for i=1:n
        x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);
    end
    %步骤2.2：令k:=k+1
    k=k+1;
    %步骤2.3：若||x_k+1-x_k||＜eps，则算法停止，输出方程组近似解x_k+1;否则令x_k:=x_k+1。
    if norm(x-x0,inf)<eps, break;end
    x0=x;
end
%步骤3：输出方程组的数值解和迭代信息。
if k>N
    Warning(['迭代次数=  ，算法超出最大迭代次数！',num2str(k)]);
else
    disp(['迭代次数=  ',num2str(k)]);
    x
end

--------------------------------------------------------------------------------------

function x = Gauss_Seidel_3_2(A,b,x0,eps,N)
%功能：用高斯-塞德尔 迭代法解n阶线性方程组Ax=b。
n=length(b);x=ones(n,1);k=0;
%步骤1：输入系数矩阵A，右端向量b,以及初始向量x0，精读eps,以及最大迭代次数N；令k=1

%步骤2：当k≤N时，执行步骤2.1~2.3。
while k<=N
    %步骤2.1：算出x1~xn的第k次迭代递推式
    for i=1:n
        x(i)=(-A(i,[1:i-1])*x([1:i-1])-A(i,[i+1:n])*x0([i+1:n])+b(i))/A(i,i);
    end
    %步骤2.2：令k:=k+1
    k=k+1;
    %步骤2.3：若||x_k+1-x_k||＜eps，则算法停止，输出方程组近似解x_k+1;否则令x_k:=x_k+1。
    if norm(x-x0,inf)<eps, break;end
    x0=x;
end
%步骤3：输出方程组的数值解和迭代信息。
if k>N
    Warning(['迭代次数=  ，算法超出最大迭代次数！',num2str(k)]);
else
    disp(['迭代次数=  ',num2str(k)]);
    x
end

5、SOR方法

%SOR.m
function x = SOR(A,b,omega,x0,eps,N)
%功能：用SOR(Successive Over-Relaxation,逐次超松弛)迭代法解n阶线性方程组Ax=b
%其中松弛因子omega(ω),初始向量x0，精度eps,最大迭代次数N由算例程序给定

%步骤1：输入矩阵A，右端向量b,松弛因子ω和初始向量x0,精度eps和最大迭代次数N。令k=0
n=length(x0);
x=ones(n,1);k=0;
%x其实可以设定为任意值，不影响后面的迭代和覆盖过程，设为全部元素为1的列向量纯粹为了表示方便

%步骤2：当k≤N时，执行步骤2.1~2.3
while k<=N
    %步骤2.1：对i=1,2,...n，计算x_i(k+1)表达式
    for i=1:n
        x(i) = omega*(-A(i,1:i-1)* x(1:i-1)-A(i,i+1:n)*x0(i+1:n)+b(i))/A(i,i)+(1-omega)*x0(i);
    end
    %步骤2.2：令k:=k+1
    k=k+1;
    %步骤2.3：若||x(k+1)-x(k)||<eps，则算法停止，输出方程组近似解x(k+1);否则令x_k:=x_k+1。
    if(norm(x-x0,inf)<eps),break;end
    x0=x;
end

%步骤3：输出信息“算法超出最大迭代次数！”，算法终止。
if k>N
    warning('算法超出最大迭代次数！');
else
    disp(['迭代次数= ',num2str(k)]);
    x
end

6、Lagrange插值法

Version 1
-----------------------------------------------------------------------------------------
%Lagrange_interpolation1.m
function b = Lagrange_interpolation1(x,y,a)
%功能：根据已知横坐标向量x和纵坐标向量y,求(n-1)次Lagrange插值多项式在点a的插值b

%步骤1：输入数据xi,yi(i=1,2,...n)及被插点a
da=a-x;%数a与向量作差，结果仍然为向量
n=length(x);%多项式的次数为n-1
l=zeros(size(y));
%步骤2：利用内置prod函数计算向量元素的乘积，并计算Lagrange(拉格朗日)插值基函数在点a的值l_j(a)
l(1)=prod(da(2:n))/prod(x(1)-x(2:n));%j=1
l(n)=prod(da(1:n-1))/prod(x(n)-x(1:n-1));%j=n
for j=2:n-1
    num=prod(da(1:j-1))*prod(da(j+1:n)); %分子
    den=prod(x(j)-x(1:j-1))*prod(x(j)-x(j+1:n)); %分母
    l(j)=num/den;
end
%步骤3：计算Lagrange插值多项式在点a处的值
b=sum(y.*l);
-----------------------------------------------------------------------------------------
version 2
-----------------------------------------------------------------------------------------
%Lagrange_interpolation2.m
function [L,p,b] = Lagrange_interpolation2(X,Y,a)
%功能：根据已知横坐标向量X和纵坐标向量Y求解n次Lagrange插值多项式L和对应的基函数p，并将a的值代入，求解对应的插值b

%步骤1：给出(n+1)个节点(xi,yi),(i=0,1,2,...,n)及插值点a，由算例程序提供
if(length(X)==length(Y))
    n=length(X);
else
    disp('X和Y的维数不相等！');
    return;
end
n=length(X);L=ones(n,n);

%步骤2：令k=1,当k≤n时，执行步骤2.1~2.3
for k=1:n
    %步骤2.1：令p_k(x):=1
    V=1;
    for i=1:n
    %步骤2.2：令i=1，对于i≤n，执行步骤2.2.1~2.2.2
        %步骤2.2.1：若i≠n，则p_k(x):=p_k(x)*(x-x_i)/(x_k-x_i)
        %步骤2.2.2：令i:=i+1
        if k~=i
            V = conv(V,poly(X(i)))/(X(k)-X(i));
        end
    end
    p(k,:)=poly2sym(V);%利用向量V创建多项式
end
p=collect(p);
p=vpa(p,4);
L=vpa(Y*p,4);
b=polyval(sym2poly(L),a);

%备注：
%1、
%（1）p = poly(r)（其中 r 是矢量）返回多项式的系数，其中多项式的根是 r 的元素。
%（2）p = poly(A)（其中 A 是 n×n 矩阵）返回矩阵 det(λI – A) 的特征多项式的 n+1 个系数
%（3）w = conv(u,v) 返回矢量 u 和 v 的卷积。如果 u 和 v 是多项式系数的矢量，对其卷积与将这两个多项式相乘等效。
%2、collect(s)函数命令使用findsym函数规定的默认变量代替collect(s,v)函数中的v;而该命令用于将符号矩阵S中所有同类项合并，并以v为符号变量输出
%3、vpa全称为variable-precision arithmetic,即算数精度。一般来说，vpa专门用来计算符号变量和函数的值。
%4、polyval(p,x)函数返回多项式p在x处的函数值，其中输入变量p=[p0 p1 p2...pn]是一个长度为n+1的横向量。

7、最小二乘法（只给出拟合函数为一次函数的程序）

function f=LeastSquare_linearfunc(X,Y,w)
%程序功能：根据已知n个点以及频数ω_i，求解最小二乘拟合曲线

%步骤1：给出n个节点(x_i,y_i)以及对应的频数w_i(算例程序给出)
%步骤2：令矩阵A为2*2的零矩阵，B为2维列零向量，令i=1,对i≤n作步骤2.1~2.2
A=zeros(2,2); B=zeros(2,1); n=length(X);
%%下面算法的解释：
%令拟合线性函数曲线为φ(x)=c1+c2*x，此时φ0(x)=1，φ1(x)=x,令法方程组为Ac=[a11,a12;a21,a22][c1;c2]=b=[b1;b2],
%则a_11=(φ0,φ0)=Σ{i=1,n}ω_iφ0(x_i)φ0(x_i)=Σ{i=1,n}ω_i
%  a_12=(φ0,φ1)=Σ{i=1,n}ω_iφ0(x_i)φ1(x_i)=Σ{i=1,n}ω_i*x_i=(φ1,φ0)=a_21
%  a_22=(φ1,φ1)=Σ{i=1,n}ω_iφ1(x_i)φ1(x_i)=Σ{i=1,n}ω_i*x_i^2
%  b_1=(f,φ0)=Σ{i=1,n}ω_i*f(x_i)*φ0(x_i)=Σ{i=1,n}ω_i*y_i
%  b_2=(f,φ1)=Σ{i=1,n}ω_i*f(x_i)*φ1(x_i)=Σ{i=1,n}ω_i*y_i*x_
for i=1:n
    %步骤2.1：算出A和B各个元素一次累加的值;%步骤2.2 令i:=i+1
    A(1,1)=A(1,1)+w(i); A(1,2)=A(1,2)+w(i)*X(i);A(2,1)=A(1,2);A(2,2)=A(2,2)+w(i)*X(i)*X(i);
    B(1,1)=B(1,1)+w(i)*Y(i);B(2,1)=B(2,1)+w(i)*Y(i)*X(i);
end
c=inv(A)*B;
syms x;
f=c(1)+c(2)*x;
f=vpa(f,4);%设置算数精度显示4位有效数字

8、不动点迭代法及其改进

%Fixed_point_Iteration_forNonlinearEquation.m
function [k,x]=Fixed_point_Iteration_forNonlinearEquation(phi,x0,eps,N)
%功能：用不动点迭代法解非线性方程φ(x)=x
%步骤1：给出初试近似值x0,容限ε，最大迭代次数N，令k=0
%步骤2：当k≤N时，执行步骤2.1~2.4
    %步骤2.1：令x=φ（x0)
    %步骤2.2：若|x-x0|＜ε,则输出x，终止
    %步骤2.3：令k:=k+1
    %步骤2.4：令x0:=x
fprintf('k    x\n');%fprintf函数：将数据写入文本文件
for k=1:N
    x=phi(x0);%phi为编程人员自定义的函数，由算例程序给出
    fprintf('%3d,  %10.9f\n',k,x) 
    %10.9f 指定输出中每行的第2个值为浮点数，字段宽度为10位数，包括小数点后的9位数

    if abs(x-x0)<eps
        return
    else
        x0=x;
        if k==N
            warning('算法超过最大迭代次数!');
        end
    end
end
----------------------------------------------------------------------------------
%Steffensen_Iteration.m
%功能：用Steffensen迭代法(不动点迭代的Aitkin加速公式改进方法)
function [k,x]=Steffenson_Iteration(phi,x0,eps,N)
fprintf('k    x\n');
for k=1:N
    y=phi(x0);
    z=phi(y);
    x=x0-(y-x0)^2/(z-2*y+x0);
    fprintf('%3d,    %10.9f\n',k,x);
    if abs(x-x0)<eps
        return
    else
        x0=x;
        if k==N
            warning('算法超过最大迭代次数!');
        end
    end
end

9、牛顿迭代法

%EquaNewton.m
function [k,x,f_value] = EquaNewton(f,Df,x0,eps,N)
%功能：用Newton(牛顿-拉夫逊)迭代法求解方程f(x)=0
fprintf('k   x   f\n');
f_value=f(x0);
for k=1:N
    x=x0-f_value/Df(x0);%Df存储f(x)的导函数
    f_value=f(x);
    fprintf( '%3d, %12.11f, %12.11f\n',k,x,f_value);
    if abs(x-x0)<eps
        return
    else
        x0=x;
        if k==N
            warning('算法超出最大迭代次数！')
        end
    end
end

10、常微分方程数值解法之改进的欧拉公式

%ImprovedEuler.m

%一阶常微分方程数值解法之改进的欧拉公式：
%步骤1：输入求解区间[a,b]，步长h，初值u0和函数f(x,u)
%步骤2：用改进的欧拉公式循环计算
    %步骤2.1：采用Euler公式提供预测值：u_(k+1)_sup(0)=u_k+hf(x_k,u_k);
    %步骤2.2：采用梯形公式计算：u_(k+1)=u_k+h/2[f(x_k,u_k)+f(x_(k+1))+u_(k+1)_sup(0)]
%步骤3：输出结果

function s=ImprovedEuler_9_1(fun,a,b,h,u0)
%[a,b]为求解区间，h为步长，s为返回向量，u0为初值
    syms x u; %定义符号变量
    n=(b-a)/h; s=zeros(1,n+1); s(1)=u0;
    for k1=1:n
        fk=double(subs(fun,{x,u},{a+(k1-1)*h,s(k1)})); %fk=f(x_k,u_k)
        fk1= double(subs(fun,{x,u},{a+k1*h,s(k1)+h*fk})); %fk1=f(x_(k+1),u_k+hf(x_k,u_k))
        s(k1+1)=s(k1)+h/2*(fk+fk1);
    end

%注：
%1、subs函数作用于第一个参数中的指定函数fun，将第三个参数（标量或向量）一一对应赋值给第二个参数中的符号变量（组）。此处即相当于每步迭代后对当前节点x_k和算得的u_k+1覆盖更新
%2、算法步骤中，下划线后面的数字代表迭代次数下标，sup(0)为上标以示区分，代表每次迭代过程中消除隐式变量u_k+1的中间变量