#include  
#include  
#include  

using namespace std;  

//求積：梯形公式  
double Func_Integral_Trapezoid  
(double lo, double hi, double(*Func)(double), int n = 1000)  
{//lo:下限；hi:上限；Func:函數(shù)；n:等分?jǐn)?shù)  
   if (n <= 0) n = 100;  
      double x;  
    double step = (hi - lo) / n;  

    double result = 0.0;  
    x = lo;  
    for (int i = 1; i < n; i++) {  
        x += step;  
        result += Func(x);  
    }  
    result += (Func(lo) + Func(hi)) / 2;  
    result *= step;  

    return result;  
}  

//求積：Simpson公式  
double Func_Integral_Simpson  
(double lo, double hi, double(*Func)(double), int n = 1000)  
{//lo:下限；hi:上限；Func:函數(shù)；n:等分?jǐn)?shù)  

    if (n <= 0) n = 100;  

    double x;  
    double step = (hi - lo) / n;  

    double result1 = 0.0;  
    x = lo;  
    for (int i = 1; i < n; i++) {  
        x += step;  
        result1 += Func(x);  
    }  
    result1 *= 2;  

    double result2 = 0.0;  
    x = lo + step / 2;  
    for (int i = 0; i < n; i++) {  
        result2 += Func(x);  
        x += step;  
    }  
    result2 *= 4;  

    double result = result1 + result2 + Func(lo) + Func(hi);  
    result *= step / 6;  

    return result;  
}  

//求積：Cotes公式  
double Func_Integral_Cotes  
(double lo, double hi, double(*Func)(double), int n = 1000)  
{//lo:下限；hi:上限；Func:函數(shù)；n:等分?jǐn)?shù)  

    if (n <= 0) n = 100;  

    double x;  
    double step = (hi - lo) / n;  

    double result1 = 0.0;  
    x = lo;  
    for (int i = 1; i < n; i++) {  
        x += step;  
        result1 += Func(x);  
    }  
    result1 *= 14;  

    double result2 = 0.0;  
    x = lo + step / 2;  

    double result3 = 0.0;  
    double x1 = lo + step / 4;  
    double x2 = lo + step / 4 * 3;  

    for (int i = 0; i < n; i++) {  
        result2 += Func(x);  
        result3 += Func(x1) + Func(x2);  
        x += step;  
        x1 += step;  
        x2 += step;  
    }  
    result2 *= 12;  
    result3 *= 32;  

    double result4 = (Func(lo) + Func(hi)) * 7;  

    double result = result1 + result2 + result3 + result4;  
    result *= step / 90;  

    return result;  
}  

//求積：Romberg公式  
double Func_Integral_Romberg  
(double lo, double hi, double(*Func)(double), int k = 4)  
{//lo:下限；hi:上限；Func:函數(shù)；k:等分指數(shù)(區(qū)間等分成2^k份)  
    int size = k + 1;  
    double *matrix = new double[size*size];  
    for (int i = 0; i < size*size; i++) matrix[i] = 0.0;  
    double step = hi - lo;  
    matrix[0] = Func_Integral_Trapezoid(lo, hi, Func, 1);  
    for (int i = 1; i < size; i++) {  
        int n = 1 << (i - 1);  
        for (int k = 0; k < n; k++) {  
            matrix[i*size + 0] += Func(lo + (k + 0.5)*step);  
        }  
        matrix[i*size + 0] *= step;  
        matrix[i*size + 0] += matrix[(i - 1)*size];  
        matrix[i*size + 0] /= 2.0;  
        step /= 2.0;  
    }  

    double temp = 1.0;  
    double factor1, factor2;  
    for (int j = 1; j < size; j++) {  
        temp *= 4.0;  
        factor1 = temp / (temp - 1);  
        factor2 = 1 / (temp - 1);  
        for (int i = j; i < size; i++) {  
            matrix[i*size + j] = factor1*matrix[i*size + j - 1]  
                - factor2*matrix[(i - 1)*size + j - 1];  
        }  
    }  

    double result = matrix[k*size + k];  
    delete[] matrix;  

    return result;  
}  

//測試用的被積函數(shù)，0到1積分為PI  
double Func_test1(double x)  
{  
    return 4 / (1 + x*x);  
}  

int main() {  
    cout << "Trapezoid Numerical Integration" << endl;  
    cout << setprecision(15) << Func_Integral_Trapezoid(0, 1, Func_test1, 1024) << endl << endl;  
    cout << "Simpson Numerical Integration" << endl;  
    cout << setprecision(15) << Func_Integral_Simpson(0, 1, Func_test1, 1024) << endl << endl;  
    cout << "Cotes Numerical Integration" << endl;  
    cout << setprecision(15) << Func_Integral_Cotes(0, 1, Func_test1, 1024) << endl << endl;  
    cout << "Romberg Numerical Integraion" << endl;  
    cout << setprecision(15) << Func_Integral_Romberg(0, 1, Func_test1, 10) << endl << endl;  
    getchar();  
    return 0;  
}

另外有需要云服務(wù)器可以了解下創(chuàng)新互聯(lián)scvps.cn，海內(nèi)外云服務(wù)器15元起步，三天無理由+7*72小時(shí)售后在線，公司持有idc許可證，提供“云服務(wù)器、裸金屬服務(wù)器、高防服務(wù)器、香港服務(wù)器、美國服務(wù)器、虛擬主機(jī)、免備案服務(wù)器”等云主機(jī)租用服務(wù)以及企業(yè)上云的綜合解決方案，具有“安全穩(wěn)定、簡單易用、服務(wù)可用性高、性價(jià)比高”等特點(diǎn)與優(yōu)勢，專為企業(yè)上云打造定制，能夠滿足用戶豐富、多元化的應(yīng)用場景需求。

創(chuàng)新互聯(lián)主要從事成都網(wǎng)站設(shè)計(jì)、成都做網(wǎng)站、網(wǎng)頁設(shè)計(jì)、企業(yè)做網(wǎng)站、公司建網(wǎng)站等業(yè)務(wù)。立足成都服務(wù)鞏留,10年網(wǎng)站建設(shè)經(jīng)驗(yàn),價(jià)格優(yōu)惠、服務(wù)專業(yè),歡迎來電咨詢建站服務(wù):13518219792
文章名稱：一重積分的C++實(shí)現(xiàn)-創(chuàng)新互聯(lián)
瀏覽地址：http://weahome.cn/article/jhcjd.html