11.2.5 Tutorial:Calculating Double Integral with NAG Functions From Origin C

1 Summary
2 What you will learn
3 Example and Steps
- 3.1 Using nag_multid_quad_monte_carlo_1
- 3.2 Using nag_multid_quad_adapt_1

Summary

Origin/OriginPro includes the complete NAG Mark 9 numerical library. This library provides multiple method for numerical integration. All functions are accessible from Origin C. In this tutorial, you will call NAG functions to perform the double integration. Note that an infinite NAG integrator is used here.

What you will learn

This tutorial will show you how to:

Call NAG functions in Origin C
Solve definite double Integral using the NAG integration routine

Example and Steps

To calculate the double integration. You need to call function in NAG library category D01 Quadrature. This category provides functions for the numerical evaluation of definite integrals in one or more dimensions. Two functions with different algorithm are available for double integral, which is nag_multid_quad_adapt_1 and nag_multid_quad_monte_carlo_1. In this tutorial, you will learn to use them to solve the integral formula below:

$\int_{0}^{1}\int_{0}^{1}xy(x+y)dxdy$

You are recommended to read the relevant tutorial Calling NAG Functions From Origin C, to study how to run and test the code below for calculation. Then you can copy and pase the code into new created .c file in Code Builder and run the compile and build process. Below is the Origin C code with comments:

Using nag_multid_quad_monte_carlo_1

#include <Origin.h>
#include <OC_nag.h> 

#define MAXCLS 20000   //maximum number of integrand evaluations to be allowed

double NAG_CALL f(Integer ndim, double x[], Nag_User *comm)  
{
  return x[0]*x[1]*(x[0]+x[1]);  //define the function formula  
}

int multid_quad_monte_carlo()
{
	Integer exit_status = 0, k, maxcls = MAXCLS, mincls;
	Integer ndim =2;  // the number of dimensions of the integral
	
	NagError fail;
	Nag_MCMethod method;
	Nag_Start cont;
	Nag_User comm;
	double a[2], b[2], acc, *comm_arr, eps, finest;		
	comm_arr=NULL;
	
	if (ndim < 1){
		printf("Invalid ndim.\n");
		exit_status = -1;
		return exit_status;   
    }	
	for (k = 0; k < ndim; k++){
		a[k] = 0.0;   // the lower limits of integration
		b[k] = 1.0;   // the upper limits of integration
	}
	
	eps = 0.01;  //the relative accuracy required
	mincls = 1000; //minimum number of integrand evaluations to be allowed
	method = Nag_ManyIterations;
	cont = Nag_Cold;
	
       /* nag_multid_quad_monte_carlo_1 (d01xbc).
	* Multi-dimensional quadrature, using Monte Carlo method,
	* thread-safe
	*/	
nag_multid_quad_monte_carlo_1(ndim, f, method, cont, a, b, &mincls, maxcls,eps, &finest, &acc, &comm_arr, &comm, &fail);

if (fail.code == NE_NOERROR || fail.code == NE_QUAD_MAX_INTEGRAND_EVAL){
	if (fail.code == NE_QUAD_MAX_INTEGRAND_EVAL){
		printf("Error from nag_multid_quad_monte_carlo_1 (d01xbc).\n%s\n",fail.message);
		exit_status = 2;
	}
	//output the calculation results
	printf("Requested accuracy = %7.2e\n", eps);
	printf("Estimated value = %7.5f\n", finest);
	printf("Estimated accuracy = %7.2e\n", acc);
	printf("Number of evaluations = %6d\n", mincls);	
}
else{
	printf("Error from nag_multid_quad_monte_carlo_1 (d01xbc).\n%s\n",fail.message);
	exit_status = 1;
}
/* Free memory allocated internally */
if (comm_arr) 
	NAG_FREE(comm_arr);

return exit_status;
}

Then you can call the functions in LabTalk Console. The results will be like this:

Requested accuracy = 1.00e-002
Estimated value = 0.33326
Estimated accuracy = 2.23e-004
Number of evaluations =   1552

Using nag_multid_quad_adapt_1

#include <OC_nag.h> 
#include <Origin.h>
#define NDIM 2   //the number of dimensions of the integral
#define MAXPTS 1000*NDIM   //maximum number of integrand evaluations to be allowed

double NAG_CALL f(Integer n, double x[], Nag_User *comm)
{
        return x[0]*x[1]*(x[0]+x[1]);  //define the function formula  
}

int  multid_quad_adapt()
{
	Integer exit_status = 0;
	Integer ndim = NDIM;
	Integer maxpts = MAXPTS;
	
	double a[2], b[2];
	Integer k;
	double finval;
	Integer minpts;
	double acc, eps;
	Nag_User comm;
	NagError fail;	
	
	for (k = 0; k < 2; k++)
	{
		a[k] = 0.0; //the lower limits of integration
		b[k] = 1.0; //the upper limits of integration
	}	
	eps = 0.001;
	minpts = 0;	

        nag_multid_quad_adapt_1(ndim, f, a, b, &minpts, maxpts, eps, &finval,&acc, &comm, &fail);

        if (fail.code != NE_NOERROR && fail.code != NE_QUAD_MAX_INTEGRAND_EVAL)
        {	
	        printf("Error from nag_multid_quad_monte_carlo_1 (d01xbc).\n%s\n",fail.message);
	        exit_status = 1;
	        return exit_status;
        }	
	//output the calculation results
	printf("Requested accuracy = %7.2e\n", eps);
	printf("Estimated value = %7.5f\n", finval);
	printf("Estimated accuracy = %7.2e\n", acc);	
	return 0;
}