Numerical Integration
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: QUAD */
/* TITLE: Numerical Integration */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: IML */
/* SUPPORT: GHG UPDATE: */
/* REF: */
/* MISC: */
/* */
/****************************************************************/
PROC IML;
/*-------------------------------------------------------------*/
/* EXAMPLE1: A SIMPLE PROGRAM THAT VERIFIES THE ADDITIVE */
/* PROPERTY OF THE INTEGRATION */
/* THE USER MUST SUPPLY */
/* 1) THE INTERVAL. */
/* 2) THE INTEGRAND . */
/* 3) OPTIONAL INTEGRATION PARAMTERS */
/* */
/* RESULTS ARE */
/* Z1 = 8.217724400583900000E-01 */
/* Z2 = 8.217724399465600000E-01 */
/*-------------------------------------------------------------*/
START FUN(T); /* DEFINE THE INTEGRAND */
V = EXP(-T*T/2.)/(1+T*T);
RETURN(V);
FINISH;
/*-------------------------------------------------------------*/
/* FIRST RUN: ON THE SEMI-INFINITE INTERVAL */
/*-------------------------------------------------------------*/
EPS = 1.E-8; /* DEFINE THE PRECISION */
A = { 0 .P }; /* DEFINE THE INTERVAL */
CALL QUAD(Z1,"FUN",A,EPS); /* CALL THE INTEGRATOR */
PRINT Z1[FORMAT =E25.15]; /* PRINT THE RESULT */
/*-------------------------------------------------------------*/
/* SECOND RUN: ON SUBINTERVALS COVERING */
/* THE SEMI-INFINITE INTERVAL */
/* WE CAN CALCULATE THE INDIVIDUAL INTEGRATIONS */
/* 1) FROM 0 TO 1 */
/* 2) FROM 1 TO 2 */
/* 3) FROM 2 TO 3 */
/* 4) FROM 3 TO INFINITY */
/* THEN WE CAN ADD THE RESULTS TO GET THE INTEGRATION */
/* FROM 0 TO INFINITY */
/*-------------------------------------------------------------*/
B = { 0 1 2 3 .P}; /* DEFINE THE INTERVAL */
CALL QUAD(Z,"FUN",B,EPS); /* CALL THE INTEGRATOR */
Z2 = SUM(Z); /* ADD THE SUBINTERVALS */
PRINT Z2[FORMAT =E25.15]; /* Z1 AND Z2 MUST BE CLOSE */
/*-------------------------------------------------------------*/
/* EXAMPLE2: A FUNCTION USEFUL FOR BIVARIATE DISTRIBUTIONS */
/* */
/* THE USER MUST SUPPLY */
/* 1) THE INTERVAL. */
/* 2) THE INTEGRAND . */
/* 3) OPTIONAL INTEGRATION PARAMTERS */
/* */
/* REF: OWEN, HANDBOOK OF STATISTICAL TABLES */
/* P 184 */
/* */
/* */
/* RESULTS ARE */
/* Z = 1.625989408808600000E-03 */
/*-------------------------------------------------------------*/
START INTEGR(Y) GLOBAL(H); /* DEFINE THE INTEGRAND */
/* PASS THE PARAMETER H */
V = EXP(-.5*H*H*(1+Y*Y))/(1+Y*Y);
RETURN(V);
FINISH;
START T_OWEN(HH,AA) GLOBAL(H,EPS);
A = AA || {.P}; /* DEFINE THE INTERVAL */
H = HH; /* POST H GLOBALLY */
CALL QUAD(Z,"INTEGR",A) EPS = EPS PEAK = AA ;
/* CALL THE INTEGRATOR */
RETURN(Z);
FINISH;
EPS = 1.E-10; /* DECLARE EPS */
Z = T_OWEN(2,1); /* EVALUATE THE FUNCTION */
PRINT Z[FORMAT =E25.15]; /* PRINT THE RESULT */
