GAMS Support Wiki

Site Tools

gams:an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule

Integrals

It is relatively easy to use some approximation scheme like Simpson's rule to calculate the value of a definite integral, One can even you 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 even use the example of an integral.