/**
变步长梯形积分    
 获取初值：  T1 = h/2［ f(x(k)) + f(x(k+1)) ］
 
                                                         n-1
 变步长梯形公式： T2n = 1/2*Tn + h/2 * ∑ f ( x(k+1/2) )
                                                        k=0
 步长： h=b-a/n

 属性： 数值积分法

《数值计算方法与算法》-2 Editon -科学出版社 P59
《C#数值计算算法编程》-周长发 P313
   
 代码维护：2007.04.20   pengkuny
**/

#include<iostream>
#include<cmath>

using namespace std;

#define f(x) (sin(x))  //举例函数
#define epsilon 0.00001  //精度

//变步长复化梯形公式
double computerAutoT(double aa, double bb)
{ 
    //迭代初值
    long n = 1;
    double h = bb-aa; //步长
    double temp = h*(f(aa) + f(bb))/2.0, t;
    double p = epsilon + 1.0;//精度控制
    double sum, x;
    
    while (p >= epsilon)
    {
        sum = 0.0;
        for (long k=0; k<n; k++)
        {
            x = aa + (k+0.5)*h;
            sum = sum + f(x);
        }

        t = (temp + h*sum)/2.0; //key step
        p = fabs(t-temp);
        temp = t;
        n = n+n;
        h = h/2.0;
    }

    cout<<"最终分点n:"<<n<<endl;
    return (t);
}


int main()
{
    double a,b;
    cout<<"变步长复化梯形积分,请输入积分范围a,b:"<<endl;
    cin>>a>>b;

    cout<<"积分结果:"<<computerAutoT(a, b)<<endl;

    system("pause");
    return 0;
}