/**
  復(fù)化Simpson積分  
                                                                  m-1                   m-1
 復(fù)化Simpson公式： Tn = (h/3) *  { f(a)/2  +  4 * ∑ f ( x_2i  + 2 * ∑ f ( x_2i+1 ) +  f(b)/2 )} 
                                                                  i=0                     i=1
 步長(zhǎng)： h=b-a/n,  n=2*m

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

《數(shù)值計(jì)算方法與算法》-2 Editon -科學(xué)出版社 P55
   
 代碼維護(hù)：2007.04.20   pengkuny
**/ 

#include<iostream>
#include<cmath>

using namespace std;

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

//復(fù)化Simpson公式
float Simpson(float aa, float bb, long n)
{//aa,bb:端點(diǎn)  n,分區(qū)數(shù)
    if (!(n%2)) //保證n為偶數(shù)
        n = n+1;
    long m = n/2;

    float sum = 0;
    float h = (bb-aa)/n; //步長(zhǎng)

    float sum1 = 0, sum2 = 0;
    for (long i=0; i<m; i++)
    {
        sum1 += f(aa + (2*i+1)*h);
    }
    for (long i=1; i<m; i++)
    {
        sum2 += f(aa + (2*i)*h);
    }

    sum = f(aa) + f(bb) + 2*sum2 + 4*sum1;

    return (h*sum/3);
}


int main()
{
    float a,b;
    long n;
    cout<<"復(fù)化Simpson積分,請(qǐng)輸入積分范圍a,b:"<<endl;
    cin>>a>>b;
    cout<<"請(qǐng)輸入分割區(qū)間數(shù)n:"<<endl;
    cin>>n;

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

    system("pause");
    return 0;
}