GSL数值积分例子

//数值计算实验 数值积分
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <gsl/gsl_integration.h>

using namespace std;

//被积函数
double f(double x)
{
    //为便于调试，先弄个有原函数的 y = x^2 + x^3 - 2*x^4
    return x*x + x*x*x - 2*x*x*x*x;
}

//被积函数 给gsl用的
double g(double x, void * params)
{
    return f(x);
}

//原函数 用于调试算法
double F(double x)
{
    return (x*x*x)/3.0 + (x*x*x*x)/4.0 - 0.4*x*x*x*x*x;
}

//输出被积函数的精确解
double Jinque(const double a, const double b)
{
    return (F(b) - F(a));
}

//梯形法 求函数在[a,b]上的定积分，积分区间分为n部分
double Tixing(const double & a, const double & b, const int & n)
{
    double sum = 0.0;
    double gaps = (b-a)/double(n);  //每个间隔的长度
    for (int i = 0; i < n; i++)
    {
        sum += (gaps/2.0) * (f(a + i*gaps) + f(a + (i+1)*gaps));
    }
    return sum;
}

//抛物线法
double Paowuxian(const double & a, const double & b, const int & n)
{
    double sum = 0.0;
    double gaps = (b-a)/double(n);  //每个间隔的长度
    double h = gaps/2.0;
    for (int i = 0; i < n; i++)
    {
        sum += (h/3.0) * (f(a + i*gaps) + f(a + (i+1)*gaps) + 4.0*f((2*a + (2*i+1)*gaps)/2.0));
    }
    return sum;
}

//柯特斯公式
double Cotes(const double & a, const double & b, const int & n)
{
    double sum = 0.0;
    double gaps = (b-a)/double(n);  //每个间隔的长度
    double h = gaps/2.0;
    for (int i = 0; i < n; i++)
    {
        sum += (h/45.0) * (7.0*f(a + i*gaps) +
                  32.0*f(a + i*gaps + 0.25*gaps) + 
                  12.0*f(a + i*gaps + 0.5*gaps) + 
                  32.0*f(a + i*gaps + 0.75*gaps) + 
                  7.0*f(a + (i+1)*gaps));
    }
    return sum;
}

//gsl解法，参考gsl文档
double gslIntegration(double & a, double & b)
{
    gsl_function gf;
    gf.function = g;

    double r, er;
    unsigned int n;
    gsl_integration_qng(&gf, a, b, 1e-10, 1e-10, &r, &er, &n);
    return r;
}

int main()
{
    double a, b;
    int n;
    cout<<"请输入积分区间:"<<endl;
    cout<<"a = ";
    cin>>a;
    cout<<"b = ";
    cin>>b;
    cout<<"请输入分割被积区间的数量:";
    cin>>n;
    if (a > b || n <= 1)
    {
        cout<<"输入错误！"<<endl;
        exit(1);
    }
    
    //设置输出精度
    cout.precision(10);

    //输出精确解
    double result = Jinque(a, b);
    cout<<"函数在["<<a<<","<<b<<"]上的定积分为:"<<result<<endl;

    //梯形法
    double result1 = Tixing(a, b, n);
    cout<<"梯形法:"<<endl;
    cout<<"函数在["<<a<<","<<b<<"]上的定积分为:"<<result1<<" 相对误差为:"
        <<abs((result1 - result)/result)*100<<"%"<<endl;

    //抛物线法
    double result2 = Paowuxian(a, b, n);
    cout<<"抛物线法:"<<endl;
    cout<<"函数在["<<a<<","<<b<<"]上的定积分为:"<<result2<<" 相对误差为:"
        <<abs((result2 - result)/result)*100<<"%"<<endl;
    
    //柯特斯公式法
    double result3 = Cotes(a, b, n);
    cout<<"柯特斯法:"<<endl;
    cout<<"函数在["<<a<<","<<b<<"]上的定积分为:"<<result3<<" 相对误差为:"
        <<abs((result3 - result)/result)*100<<"%"<<endl;
    
    //调用gsl函数
    double result4 = gslIntegration(a, b);
    cout<<"gsl函数结果:"<<endl;
    cout<<"函数在["<<a<<","<<b<<"]上的定积分为:"<<result4<<" 相对误差为:"
        <<abs((result4 - result)/result)*100<<"%"<<endl;
    return 0;
}

http://www.cnblogs.com/make217/archive/2013/04/02/2995108.html