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;
}