/**
變步長Simpson積分    
 1.獲取初值：  T1 = h/2［ f(a) + f(b) ］,n=1, 步長： h=b-a/n, 且令Sn = Tn
 
                                                                  n-1
 2.用變步長梯形公式計算： T2n = 1/2*Tn + h/2 * ∑ f ( x(k+1/2) ) 
                                                                  k=0
 3.用Simpson求積:S_2n = (4T_2n - T_n ) /3
 不滿足精度,則加倍分點n,迭代求值.

 屬性： 數值積分法

《數值計算方法與算法》-2 Editon -科學出版社 P59
《C#數值計算算法編程》-周長發 P315
   
 代碼維護：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 t1 = h*(f(aa) + f(bb))/2.0, t2;//t1表示Tn, t2表示T2n
    double s1=t1, s2=0;          //s1表示Sn, s2表示S2n
    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);
        }

        t2 = (t1 + h*sum)/2.0; //key step
        s2 = (4.0*t2 - t1)/3.0; //key step
        p = fabs(s2-s1);
        t1 = t2; s1 = s2; n = n+n; h = h/2.0;
    }

    cout<<"最終分點n:"<<n<<endl;
    return (s2);
}


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

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

    system("pause");
    return 0;
}