2-3、迭代法计算定积分

内容来自王晓华老师

这块内容有点硬核，先做了解，主要学习如何使用迭代解决问题的步骤

两种计算数值积分的常用算法，分别是变步长梯形公式法和变步长辛普森公式法。首先从梯形公式入手来推导出复合梯形公式法，在实现复合梯形公式法的基础上，再实现变步长梯形公式法。

同样，变步长辛普森公式法也是从辛普森公式入手的，首先实现复合辛普森公式法的算法，然后再实现变步长辛普森公式法

梯形公式法

复合梯形公式法
用梯形公式计算定积分，当区间 [a,b] 比较大的时候，其误差也会大到无法接受。如果将大的区间分割成 n 个小的区间，对每个小区间应用梯形公式计算内接梯形面积，然后再将所有小梯形的面积求和得到结果，误差则会大大的缩小，这就是复合梯形公式法。这样分割以后，每个小区间的长度被称为“步长”，即 step=(b-a)/n。积分区间分割得越细，梯形公式法得到的近似值就越逼近真实值。

double trapezium(std::function<double(double)> func, double a, double b, int n)
{
    double step = (b - a) / n;
    double sum = (func(a) + func(b)) / 2.0;

    for (int j = 1; j < n; j++)
    {
        double xj = a + j*step;
        sum = sum + func(xj);
    }

    sum *= step;

    return sum;
}

变步长梯形公式法
为了提高迭代计算的效率，人们想出了一种利用迭代思想的变步长梯形公式法。在对定积分图形的内接梯形分割的时候，每个迭代都把上一个迭代分割的梯形再平均分割成两个小梯形。随着迭代次数增加，逐步增加梯形分割的个数，使得梯形分割沿积分自变量方向的步长由粗到细，逐步变化，就是所谓的变步长。每个迭代计算一次小梯形面积的和，并与上个迭代计算的小梯形面积的和做比较，若相邻两次迭代的差值达到精度要求，则退出迭代计算，否则就对当前的所有小梯形继续分割，进行下次迭代计算。如图（2）所示，每次分割得到的两个小梯形的面积之和都比分割前的大梯形面积更接近曲线的积分值。

double vs_trapezium(std::function<double(double)> func, double a, double b, double eps)
{
    int n = 1;   //初始分割一个大梯形
    double t1 = (b - a) * (func(a) + func(b)) / 2.0; //用梯形公式计算初始梯形的面积
    double diff = eps + 1.0;
    do
    {
        n = 2 * n;    //梯形分割加倍
        double t = trapezium(func, a, b, n);  //用复合梯形公式法计算 n 个小梯形的面积之和
        diff = std::fabs(t1 - t);      //计算两次迭代的结果差值
        t1 = t; //更新迭代变量
    } while (diff >= eps);

    return t1;
}

辛普森公式法

复合辛普森公式法

double simpson(std::function<double(double)> func, double a, double b, int n)
{
    double s1, s2;
    int i;

    double step = (b - a) / n;
    s1 = s2 = 0;
    for (i = 1; i < n; i++)
    {
        s2 += func(a + step * i);   //xi 求和
    }
    for (i = 1; i <= n; i++)
    {
        s1 += func(a + step * i - step / 2);  //(xi - step/2)求和
    }

    double s = step * (func(a) + func(b) + 4 * s1 + 2 * s2) / 6.0;

    return s;
}

可变长辛普森公式法
和梯形公式一样，复合辛普森公式也可以改造为变步长辛普森公式法。改造的方法就是使用迭代法的思想，通过改变区间个数 n 使得步长 step 也跟着变化，当迭代差值符合精度要求时即可停止迭代。算法的迭代变量仍然是每次分割后的小区间上使用辛普森公式计算的插值曲线面积之和，迭代关系则非常简单，就是用本迭代的迭代变量代替上个迭代的迭代自变量的值，迭代终止条件就是两个迭代的迭代变量之差小于精度值。迭代变量的初始值就是在区间 [a,b] 上应用辛普森公式计算最大的区间面积。用一个变量 n 表示当前迭代分割小梯形的个数，n 的值每个迭代增加一倍。而每次分割后的小区间面积和的计算可由第 2-2 课中给出的复合辛普森公式算法 simpson() 函数计算，迭代算法的整体结构与变步长梯形法类似。

double vs_simpson(std::function<double(double)> func, double a, double b, double eps)
{
    int n = 1;   //初始分割一个大梯形
    double s1 = (b - a) * (func(a) + func(b) + 4 * func((a + b) / 2.0)) / 6.0; //计算初始梯形的面积
    double diff = eps + 1.0;
    do
    {
        n = 2 * n;    //梯形分割加倍
        double t = simpson(func, a, b, n); //用复合辛普森公式计算 n 个小梯形的面积之和
        diff = std::fabs(s1 - t);      //计算两次迭代的结果差值
        s1 = t; //更新迭代变量
    } while (diff >= eps);

    return s1;
}