/**
变步长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;
}