3.9.5.1 OriginC: Quadrature Integral
Using Origin C, user can access NAG integral routines to perform integration.
And this page will show you how to access the powerful NAG integrators to do integration on different kinds of integrands. For example, user can do integration on normal integrand, integrand with parameters, integrand with oscillation, infinite integral, higher-dimension integration etc.
Simple Integral Function(Origin80 SR4)
The following integral routine shows how to do a basic integration on a simple integrand.
The integrand has only one integration variable.
#include <OC_nag.h>
// Start your functions here.
//NAG_CALL denotes proper calling convention. You may treat it like a function pointer
//and define your own integrand
static double NAG_CALL f_callback_ex(double x, Nag_User *comm)
{
int *use_comm = (int *)comm->p;
return (x*sin(x*30.0)/sqrt(1.0-x*x/(PI*PI*4.0)));
}
void nag_d01sjc_ex()
{
double a = 0.0;
double b = PI * 2.0; // integration interval
double epsabs, abserr, epsrel, result;
// you may use epsabs and epsrel and this quantity to enhance your desired precision
// when not enough precision encountered
epsabs = 0.0;
epsrel = 0.0001;
// The max number of sub-intervals needed to evaluate the function in the integral
// The more diffcult the integrand the larger max_num_subint should be
// For most problems 200 to 500 is adequate and recommmended
int max_num_subint = 200;
Nag_QuadProgress qp;
NagError fail;
Nag_User comm;
static int use_comm[1] = {1};
comm.p = (Pointer)&use_comm;
//Note:nag_1d_quad_gen (d01ajc) has been replaced by nag_1d_quad_gen_1 (d01sjc) at Mark 25
d01sjc(f_callback_ex, a, b, epsabs, epsrel, max_num_subint,&result, &abserr, &qp, &comm, &fail);
// For the error other than the following three errors which are due to bad input parameters
// or allocation failure NE_INT_ARG_LT NE_BAD_PARAM NE_ALLOC_FAIL
// You will need to free the memory allocation before calling the integration routine again to
// avoid memory leakage
if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
{
NAG_FREE(qp.sub_int_beg_pts);
NAG_FREE(qp.sub_int_end_pts);
NAG_FREE(qp.sub_int_result);
NAG_FREE(qp.sub_int_error);
}
printf("%10.6f", result);
}
Integral Function with Parameters(Origin80 SR4)
The following example shows how to define and perform integration on a integrand with parameters.
Note that the paramters are passed to the integrator by a user struct, which can avoid using static variables as parameters of the integrand. And it is thread-safe by this way.
The example code can also be adapted to using infinite integrator of NAG. For instance, by uncomment the line using infinite integrator d01smc, the example code can then be used to perform infinite integration.
#include <OC_nag.h>
// Start your functions here.
struct user // parameters in the integrand
{
double A, Xc, W;
};
// Function supplied by user, return the value of the integrand at a given x.
static double NAG_CALL f_callback(double x, Nag_User *comm)
{
struct user *sp = (struct user *)(comm->p);
double A, xc, w; // temp variable to accept the parameters in the Nag_User communication struct
A= sp->A;
xc= sp->Xc;
w= sp->W;
return A* exp( -2*(x - xc)*(x - xc)/w/w) / (w*sqrt(PI/2));
}
void nag_d01sjc_ex()
{
double a = 0.0;
double b = 2.0; // integration interval
// Through the absolute accuracy epsabs, relative accuracy epsrel and max_num_subint you can
// control the precision of the integration you need
// if epsrel is set negative, the absolute accuracy will be used.
// Similarly, you can control only relative accuracy by set the epsabs negative
double epsabs = 0.0, epsrel = 0.0001;
// The max number of sub-intervals needed to evaluate the function in the integral
// The more diffcult the integrand the larger max_num_subint should be
// For most problems 200 to 500 is adequate and recommmended
int max_num_subint = 200;
// Result keeps the approximate integral value returned by the algorithm
// abserr is an estimate of the error which should be an upper bound for the |I - result|
// where I is the integral value
double result, abserr;
// The structure of type Nag_QuadProgress,
// it contains pointers allocated memory internally with max_num_subint elements
Nag_QuadProgress qp;
// The NAG error parameter (structure)
static NagError fail;
// Parameters passed to integrand by Nag_User communication struct
Nag_User comm;
struct user s;
double A = 1.0;
double xc = 0.0;
double w =1.0;
s.A = A;
s.Xc = xc;
s.W = w;
comm.p = (Pointer)&s;
// Perform integration
// There are 3 kinds of infinite boundary types you can use in Nag infinite integrator
// Nag_LowerSemiInfinite, Nag_UpperSemiInfinite, Nag_Infinite
//d01smc(f_callback, Nag_LowerSemiInfinite, b, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
d01sjc(f_callback, a, b, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
// you may want to exam the error by printing out error message, just uncomment the following lines
// if (fail.code != NE_NOERROR)
// printf("%s\n", fail.message);
// For the error other than the following three errors which are due to bad input parameters
// or allocation failure NE_INT_ARG_LT NE_BAD_PARAM NE_ALLOC_FAIL
// You will need to free the memory allocation before calling the integration routine again to avoid memory leakage
if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
{
NAG_FREE(qp.sub_int_beg_pts);
NAG_FREE(qp.sub_int_end_pts);
NAG_FREE(qp.sub_int_result);
NAG_FREE(qp.sub_int_error);
}
printf("%10.6f", result);
}
Multi-dimension Integral Function(Origin80 SR0)
For higher dimension (more than 2D) integral, user can access NAG integrator d01wcc to perform higher dimension integration.
#include <OC_nag.h>
#define NDIM 4
#define MAXPTS 1000*NDIM
#define EXIT_SUCCESS 1
#define EXIT_FAILURE 0
static double NAG_CALL f(Integer n, double z[], Nag_User *comm)
{
double tmp_pwr;
tmp_pwr = z[1]+1.0+z[3];
return z[0]*4.0*z[2]*z[2]*exp(z[0]*2.0*z[2])/(tmp_pwr*tmp_pwr);
}
int nag_d01wcc_ex()
{
Integer ndim = NDIM; // the integral dimension
Integer maxpts = MAXPTS; // maximum number of function evaluation
double a[4], b[4];
Integer k;
static NagError fail;
double finval;
Integer minpts;
double acc, eps;
Nag_User comm;
for (k=0; k < 4; ++k) // integration interval
{
a[k] = 0.0;
b[k] = 1.0;
}
eps = 0.0001; // set the precision
minpts = 0;
d01wcc(ndim, f, a, b, &minpts, maxpts, eps, &finval, &acc, &comm, &fail);
if (fail.code != NE_NOERROR)
printf("%s\n",fail.message);
if (fail.code == NE_NOERROR || fail.code == NE_QUAD_MAX_INTEGRAND_EVAL)
{
printf("Requested accuracy =%12.2e\n", eps);
printf("Estimated value =%12.4f\n", finval); // the final integration result
printf("Estimated accuracy =%12.2e\n", acc);
return EXIT_SUCCESS;
}
else
return EXIT_FAILURE;
}
|