Integrals
Posted: Fri Jun 05, 2020 8:20 am
How do I calculate the value of a definite integral in GAMS?
Code: Select all
* 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;