/**
變步長(zhǎng)梯形積分    
 獲取初值：  T1 = h/2［ f(x(k)) + f(x(k+1)) ］
 
                                                         n-1
 變步長(zhǎng)梯形公式： T2n = 1/2*Tn + h/2 * ∑ f ( x(k+1/2) )
                                                        k=0
 步長(zhǎng)： h=b-a/n

 屬性： 數(shù)值積分法

《數(shù)值計(jì)算方法與算法》-2 Editon -科學(xué)出版社 P59
《C#數(shù)值計(jì)算算法編程》-周長(zhǎng)發(fā) P313
   
 代碼維護(hù)：2007.04.20   pengkuny
**/

#include<iostream>
#include<cmath>

using namespace std;

#define f(x) (sin(x))  //舉例函數(shù)
#define epsilon 0.00001  //精度

//變步長(zhǎng)復(fù)化梯形公式
double computerAutoT(double aa, double bb)
{ 
    //迭代初值
    long n = 1;
    double h = bb-aa; //步長(zhǎng)
    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<<"最終分點(diǎn)n:"<<n<<endl;
    return (t);
}


int main()
{
    double a,b;
    cout<<"變步長(zhǎng)復(fù)化梯形積分,請(qǐng)輸入積分范圍a,b:"<<endl;
    cin>>a>>b;

    cout<<"積分結(jié)果:"<<computerAutoT(a, b)<<endl;

    system("pause");
    return 0;
}