test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

    int status = 0; size_t neval = 0 ;
    double result = 0, abserr = 0 ;

    double exp_result = 3.222948711817264211E+01;
    double exp_abserr = 2.782360287710622870E+01;
    int exp_neval  =  87;
    int exp_ier    =  GSL_ETOL;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-3,
                                  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qng(f1) sing beyond 87pt result");
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) sing beyond 87pt abserr");
    gsl_test_int((int)neval,exp_neval,"qng(f1) sing beyond 87pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) sing beyond 87pt status") ;

    status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-3,
                                  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse beyond 87pt result");
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) rev beyond 87pt abserr");
    gsl_test_int((int)neval,exp_neval,"qng(f1) rev beyond 87pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) rev beyond 87pt status") ;
  }

  /* Test the adaptive integrator QAG */

  {
    int status = 0, i; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    double exp_result = 7.716049382715854665E-02 ;
    double exp_abserr = 6.679384885865053037E-12 ;
    int exp_neval  =     165;
    int exp_ier    =       0;
    int exp_last   =       6;

    double a[6] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125 } ;
    double b[6] = { 0.03125, 1, 0.5, 0.25, 0.125, 0.0625 } ;
    double r[6] = { 3.966769831709074375E-06, 5.491842501998222409E-02,
                    1.909827770934243926E-02, 2.776531175604360531E-03,
                    3.280661030752063693E-04, 3.522704932261797744E-05 } ;
    double e[6] = { 6.678528276336181873E-12, 6.097169993333454062E-16,
                    2.120334764359736934E-16, 3.082568839745514608E-17,
                    3.642265412331439511E-18, 3.910988124757650942E-19 } ;
    int order[6] = { 1, 2, 3, 4, 5, 6 } ;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha) ;

    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
                                  GSL_INTEG_GAUSS15, w,
                                  &result, &abserr) ;

    gsl_test_rel(result,exp_result,1e-15,"qag(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) smooth abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f1) smooth neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f1) smooth last") ;  
    gsl_test_int(status,exp_ier,"qag(f1) smooth status") ;

    for (i = 0; i < 6 ; i++) 
        gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1) smooth alist") ;

    for (i = 0; i < 6 ; i++) 
        gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1) smooth blist") ;

    for (i = 0; i < 6 ; i++) 
        gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1) smooth rlist") ;

    for (i = 0; i < 6 ; i++) 
        gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1) smooth elist") ;

    for (i = 0; i < 6 ; i++) 
        gsl_test_int((int)w->order[i],order[i]-1,"qag(f1) smooth order") ;

    p.neval = 0;

    status = gsl_integration_qag (&fc, 1.0, 0.0, 0.0, 1e-10, w->limit,
                                  GSL_INTEG_GAUSS15, w,
                                  &result, &abserr) ;

    gsl_test_rel(result,-exp_result,1e-15,"qag(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) reverse abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f1) reverse neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f1) reverse last") ;  
    gsl_test_int(status,exp_ier,"qag(f1) reverse status") ;

    gsl_integration_workspace_free (w) ;

  }

  /* Test the same function using an absolute error bound and the
     21-point rule */

  {
    int status = 0, i; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    double exp_result = 7.716049382716050342E-02 ;
    double exp_abserr = 2.227969521869139532E-15 ;
    int exp_neval  =     315;
    int exp_ier    =       0;
    int exp_last   =       8;

    double a[8] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625,
                    0.0078125 } ;
    double b[8] = { 0.0078125, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125,
                    0.015625 } ;
    double r[8] = { 3.696942726831556522E-08, 5.491842501998223103E-02,
                    1.909827770934243579E-02, 2.776531175604360097E-03,
                    3.280661030752062609E-04, 3.522704932261797744E-05,
                    3.579060884684503576E-06, 3.507395216921808047E-07 } ;
    double e[8] = { 1.371316364034059572E-15, 6.097169993333454062E-16,
                    2.120334764359736441E-16, 3.082568839745514608E-17,
                    3.642265412331439511E-18, 3.910988124757650460E-19,
                    3.973555800712018091E-20, 3.893990926286736620E-21 } ;
    int order[8] = { 1, 2, 3, 4, 5, 6, 7, 8 } ;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, 0.0, 1.0, 1e-14, 0.0, w->limit,
                                  GSL_INTEG_GAUSS21, w,
                                  &result, &abserr) ;

    gsl_test_rel(result,exp_result,1e-15,"qag(f1,21pt) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) smooth abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) smooth neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) smooth last") ;  
    gsl_test_int(status,exp_ier,"qag(f1,21pt) smooth status") ;

    for (i = 0; i < 8 ; i++) 
        gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1,21pt) smooth alist") ;

    for (i = 0; i < 8 ; i++) 
        gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1,21pt) smooth blist") ;

    for (i = 0; i < 8 ; i++) 
        gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1,21pt) smooth rlist") ;

    for (i = 0; i < 8 ; i++) 
        gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1,21pt) smooth elist") ;

    for (i = 0; i < 8 ; i++) 
        gsl_test_int((int)w->order[i],order[i]-1,"qag(f1,21pt) smooth order");


    p.neval = 0;
    status = gsl_integration_qag (&fc, 1.0, 0.0, 1e-14, 0.0, w->limit,
                                  GSL_INTEG_GAUSS21, w,
                                  &result, &abserr) ;

    gsl_test_rel(result,-exp_result,1e-15,"qag(f1,21pt) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) reverse abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) reverse neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) reverse last") ;  
    gsl_test_int(status,exp_ier,"qag(f1,21pt) reverse status") ;

    gsl_integration_workspace_free (w) ;

  }

  /* Adaptive integration of an oscillatory function which terminates because
     of roundoff error, uses the 31-pt rule */

  {
    int status = 0; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    double exp_result = -7.238969575482959717E-01;
    double exp_abserr =  1.285805464427459261E-14;
    int exp_neval   =     31;
    int exp_ier     =     GSL_EROUND;
    int exp_last    =     1;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, 0.3, 2.71, 1e-14, 0.0, w->limit, 
                                  GSL_INTEG_GAUSS31, w, 
                                  &result, &abserr) ;

    gsl_test_rel(result,exp_result,1e-15,"qag(f3,31pt) oscill result");
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) oscill abserr");
    gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) oscill neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) oscill last") ;  
    gsl_test_int(status,exp_ier,"qag(f3,31pt) oscill status") ;

    p.neval = 0;
    status = gsl_integration_qag (&fc, 2.71, 0.3, 1e-14, 0.0, w->limit, 
                                  GSL_INTEG_GAUSS31, w, 
                                  &result, &abserr) ;

    gsl_test_rel(result,-exp_result,1e-15,"qag(f3,31pt) reverse result");
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) reverse abserr");
    gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) reverse neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) reverse last") ;  
    gsl_test_int(status,exp_ier,"qag(f3,31pt) reverse status") ;

    gsl_integration_workspace_free (w) ;

  }

  /* Check the singularity detection (singularity at x=-0.1 in this example) */

  {
    int status = 0; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    int exp_neval  =     5151;
    int exp_ier    =     GSL_ESING;
    int exp_last   =     51;

    double alpha = 2.0 ;
    gsl_function f = make_function(&f16, &alpha);

    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit,
                                  GSL_INTEG_GAUSS51, w, 
                                  &result, &abserr) ;

    gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) sing neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) sing last") ;  
    gsl_test_int(status,exp_ier,"qag(f16,51pt) sing status") ;

    p.neval = 0;
    status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit,
                                  GSL_INTEG_GAUSS51, w, 
                                  &result, &abserr) ;

    gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) rev neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) rev last") ;  
    gsl_test_int(status,exp_ier,"qag(f16,51pt) rev status") ;

    gsl_integration_workspace_free (w) ;

  }

  /* Check for hitting the iteration limit */

  {
    int status = 0, i; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (3) ;

    double exp_result =  9.565151449233894709 ;
    double exp_abserr =  1.570369823891028460E+01;
    int exp_neval  =     305;
    int exp_ier    =     GSL_EMAXITER;
    int exp_last   =     3;

    double a[3] = { -5.000000000000000000E-01,
                    0.000000000000000000,
                    -1.000000000000000000 } ;
    
    double b[3] = { 0.000000000000000000,
                    1.000000000000000000,
                    -5.000000000000000000E-01 } ;
    
    double r[3] = { 9.460353469435913709,
                    9.090909090909091161E-02,
                    1.388888888888888812E-02 } ;
    
    double e[3] = { 1.570369823891028460E+01,
                    1.009293658750142399E-15,
                    1.541976423090495140E-16 } ;
    
    int order[3] = { 1, 2, 3 } ;

    double alpha = 1.0 ;
    gsl_function f = make_function(&f16, &alpha);
    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit, 
                                  GSL_INTEG_GAUSS61, w, 
                                  &result, &abserr) ;

    gsl_test_rel(result,exp_result,1e-15,"qag(f16,61pt) limit result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) limit abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) limit neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) limit last") ;  
    gsl_test_int(status,exp_ier,"qag(f16,61pt) limit status") ;

    for (i = 0; i < 3 ; i++) 
        gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f16,61pt) limit alist") ;

    for (i = 0; i < 3 ; i++) 
        gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f16,61pt) limit blist") ;

    for (i = 0; i < 3 ; i++) 
        gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f16,61pt) limit rlist") ;

    for (i = 0; i < 3 ; i++) 
        gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f16,61pt) limit elist") ;

    for (i = 0; i < 3 ; i++) 
        gsl_test_int((int)w->order[i],order[i]-1,"qag(f16,61pt) limit order");

    p.neval = 0;
    status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit, 
                                  GSL_INTEG_GAUSS61, w, 
                                  &result, &abserr) ;

    gsl_test_rel(result,-exp_result,1e-15,"qag(f16,61pt) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) reverse abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) reverse neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) reverse last") ;  
    gsl_test_int(status,exp_ier,"qag(f16,61pt) reverse status") ;