## 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        */

```