7.5　数值积分

STL为很多复杂的算法提供基本的计算。使用STL，程序员可以方便地实现这些算法。下面以数值积分为例介绍STL对算法的支持。积分的思想是使用生成器产生一组点序列，其中生成器（generator）是一个类，它重载了函数调用运算符operator()。STL中的算法

generate(iterator b, iterator e, generator g)

为函数产生(0,1)区间的vector的值。程序中要用到algorithm、numeric和vector库。

文件 stl_integration1.cpp

// Simple integration routine for x * x over (0, 1)

#include <iostream>

#include <numeric>

#include <algorithm>

#include <vector>

using namespace std;

// The function is represented in class gen

class gen {             // generator for integrated function

public:

   gen(double x_zero, double increment)

      : x(x_zero), incr(increment) { }

      double operator()() { x += incr; return x * x; }

private:

   double x, incr;

};

double integrate(gen g, int n) // integrate on (0,1)

{

   vector<double> fx(n);

   generate(fx.begin(), fx.end(), g);

   return(accumulate(fx.begin(), fx.end(), 0.0) / n);

}

int main()

{

   const int n = 10000;

   gen g(0.0, 1.0/n);

   cout << "integration program x**2" << endl;

   cout << integrate(g, n) << endl;

}

stl_integration 程序解析

■ class gen { // generator for integrated function

 public:

    gen(double x_zero, double increment)

    : x(x_zero), incr(increment) { }

    double operator()() {x += incr; return x * x; }

 private:

    double x, incr;

 };

要编写一个生成器必须编写一个类来重载无参数的函数调用运算符。本例中gen::operator()()用来增加x的值，类构造函数用来初始化变量x以及增量incr。这里的基本思想是使用生成器生成一个值序列，使用这些值进行数值积分。

■ double integrate(gen g, int n) // integrate (0,1)

 {

     vector<double> fx(n);

我们使用生成器在0~1以1/n为步长进行积分。表达式x*x的值存储在向量fx中，其值在0~1（包括1）。

■ generate(fx.begin(), fx.end(), g);

生成器对象g用于以1/n为步长生成n x*x的值。

■ return(accumulate(fx.begin(), fx.end(), 0.0) / n);

这个数值累加和是积分f(x) = x2的近似值。

■ int main()

 {

    const int n = 10000;

    gen g(0.0, 1.0/n);

    cout << "integration program x**2" << endl;

    cout << integrate(g, n) << endl;

 }

生成器对象的初始步长为1/10 000，这意味着求值时要使用10 000个点。这个程序在我们的系统中的答案为0.333 383，而使用微积分计算得到的答案是1/3，因此数值解法可以精确到小数点后4位。

可以编写一个更精确的数值积分器。我们使用一系列矩形近似求取曲线围成的面积，其中，增量是矩形的宽，函数值是矩形的高。计算矩形高度时，增量提供了两个选择。前面的例子中使用矩形右边的点计算高度，而这个例子中为了提高积分的数值精度，分别使用矩形两边的点计算最小高度和最大高度，然后再分别使用这两个高度计算矩形的面积，最后取面积的平均值。

文件 stl_integration2.cpp

// Simple integration routine for x * x over (0, 1)

#include <iostream>

#include <numeric>

#include <algorithm>

#include <vector>

using namespace std;

// The function is represented in class gen

class gen {     // generator for integrated function

public:

   gen(double x_zero, double increment)

    : x(x_zero), incr(increment) { }

   double operator()() { x += incr; return x * x; }

private:

   double x, incr;

};

// Integrate on (0,1)

double integrate(gen g, int n, double& diff)

{

   vector<double> fx(n + 1), sm(n), lg(n);

   double s, l;

   generate(fx.begin(), fx.end(), g);

   for (int i = 0; i < n; ++i)

   if (fx[i] > fx[i + 1]) {

      sm[i] = fx[i + 1]; lg[i] = fx[i];

   }

   else {

      sm[i] = fx[i]; lg[i] = fx[i + 1];

   }

   s = accumulate(sm.begin(), sm.end(), 0.0) / n;

   l = accumulate(lg.begin(), lg.end(), 0.0) / n;

   diff = l - s;

   return (s + l) / 2;

}

int main()

{

   const int n = 10000;

   gen g(0.0, 1.0/n);

   cout << "integration program x**2" << endl;

   double d, i = integrate(g, n, d);

   cout << "integral = " << i << " +/- "

        << (d / 2) << endl;

}

第二个版本的integrate()的计算更精确，它对diff进行了误差估算。这个方法在我们的系统上运行的结果是0.333 333，可以精确到小数点后6位。