一重積分的C++實現

#include<iostream>  
#include<iomanip>  
#include<cmath>  

using namespace std;  

//求積：梯形公式  
double Func_Integral_Trapezoid  
(double lo, double hi, double(*Func)(double), int n = 1000)  
{//lo:下限；hi:上限；Func:函數；n:等分數  
   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:函數；n:等分數  

    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:函數；n:等分數  

    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:函數；k:等分指數(區間等分成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;  
}  

//測試用的被積函數，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;  
}