How do I calculate the value of a definite integral in GAMS?

It is relatively easy to use some approximation scheme like Simpson’s rule to calculate the value of a definite integral. One can use the GAMS macro facility to make this look decent:

```
* Composite Simpson rule
$macro simpson(f,a,b,nn) \
abort$(mod(nn,2) or nn<>round(nn)) 'nn not an even integer'; \
dx = (b-a)/nn; \
x = a + dx; \
interval = f(a) + f(b); \
m = 4; \
for(ii = 2 to nn, \
interval = interval + m*f(x); \
m = 6 - m; \
x = x + dx ); \
interval = dx*interval/3;
* following will be used by the macro simpson
scalar n 'number of intervals - has to be even'
a 'interval low'
b 'interval high'
dx 'sub interval'
m 'constant'
ii 'loop index'
x 'function argument'
interval 'result';
* now we use a known integral to test the macro
* integral of sin(x) for [0,1] = 1 - cos(1);
scalar exactResult, diff; exactResult = 1-cos(1); display exactResult;
simpson(sin,0,1,2); diff = abs(interval-exactResult); display diff;
simpson(sin,0,1,10); diff = abs(interval-exactResult); display diff;
simpson(sin,0,1,1000); diff = abs(interval-exactResult); display diff;
* now we define a macro for the function
$macro g(y) sin(y)
simpson(g,0,1,10); diff = abs(interval-exactResult); display diff;
* and an even more complicated one - normal density
$macro p(x) 1/(sigma*sqrt(2*PI))*exp(-sqr(x-mu)/(2*sqr(sigma)))
parameters sigma, mu;
set i / 1*20 /,
j / 1*25 /;
parameters sx(i), mx(i,j), result(i,j);
sx(i) = uniform(1,2); mx(i,j) = uniform(-1,+1);
loop((i,j),
sigma = sx(i); mu = mx(i,j);
simpson(p,0,10,100);
result(i,j) = interval;
);
display result;
```

This implementation has the disadvantage that it can’t be used in model equations. GAMS has the ability to use extrinsic functions and external equations to implement some integral function that can be even used in model algebra. The discussion about the differences between the two methods includes the example of an integral.