演算法基本原理：把原區間分為一系列小期間(n份)，在每個小區間上都用小的梯形面積來近似代替原函式的積分，當小區間足夠小時，就可以得到原來積分的近似值。但是這個過程中有一個問題是步長的取值，步長太大精度難以保證，步長太小會導致計算量的增加，所以，實際計算中常常採用變步長的方法，使得在步長逐次減小的過程中，求得的積分結果滿足要求的精度為止。

       首先，給出兩個計算公式

              （1）

              （2） //將步長h二分，計算以h/2為步長的積分

步驟：

• 取n=1，利用公式（1）計算積分值
• 進行二分，利用新的公式（2）計算新的積分值
• 進行判斷，如果兩次計算積分值的差在給定的誤差範圍之內，二分後的積分值就是所需的結果，計算結束，否則返回第二步繼續執行

因此，主要計算的問題有兩個

一. 被積函式值的運算（設定Function類）

二. 變步長梯形積分的實現（設定Trapz類）

#ifndef _TRAPZINT_H_
#define _TRAPZINT_H_

class Function{
public:
	virtual double operator()(double x) const = 0;//純虛擬函式過載運算子()
	virtual ~Function(){}//解構函式
};

class MyFunction :public Function{
public:
	virtual double operator()(double x) const;//覆蓋虛擬函式
};

class Integration{//抽象類Integration的定義
public:
	virtual double operator()(double a, double b, double eps) const = 0;
	virtual ~Integration(){}
};

class Trapz :public Integration{//公有派生類Trapz的定義
public:
	Trapz(const Function &f) :f(f){}//建構函式
	virtual double operator()(double a, double b, double eps) const;//覆蓋虛擬函式
private:
	const Function &f;//私有成員，Function類物件的指標
};

#endif
 

#include"Trapzint.h"
#include<iostream>
#include<cmath>
using namespace std;

//計算被積函式的函式值
double MyFunction::operator()(double x) const{
	return log(1.0 + x) / (1.0 + x*x);
}

//積分運算過程，過載為運算子()
double Trapz::operator()(double a, double b, double eps) const{
	//計算區間[a,b]之間誤差在eps之內的積分值
	bool done = false;//Trapz類的虛擬函式成員
	int n = 1;
	double h = b - a;
	double tn = h*(f(a) + f(b)) / 2;//計算n=1時的積分
	double t2n;//步長為一半的時候的積分值
	do{
		double sum = 0;
		for (int k = 0; k < n; k++){
			double x = a + (k + 0.5)*h;
			sum += f(x);
		}
		t2n = (tn + h*sum) / 2.0;//利用公式（2）計算積分
		if (fabs(tn - t2n) < eps)
			done = true;//判斷積分誤差
		else{
			tn = t2n;
			n *= 2;
			h /= 2;
		}
	} while (!done);
	return t2n;
}

#include"Trapzint.h"
#include<iostream>
#include<iomanip>
using namespace std;
int main(){
	MyFunction f;//定義MyFunction類的物件
	Trapz trapz(f);//定義Trapz類的物件
	//計算並輸出結果
	cout << "Trapz Int:" << setprecision(7) << trapz(0, 2, 1e-7) << endl;
	return 0;
}

結果