User Tools

Site Tools


gams:an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule

An Example how to use the GAMS $macro feature to implement the Simpson's rule

* Composit Simpson rule - for details see Simpson's rule at Wikipedia (http://en.wikipedia.org/wiki/Simpson%27s_rule)
$macro simpson(f,a,b,nn)     \
       n  = max(2,ceil(nn)); \
       n = n + mod(n,2);     \
       dx = (b-a)/n;         \
       x = a + dx;           \
       int = f(a) + f(b);    \
       m = 4;                \
       for(ii = 2 to n,      \
         int = int + m*f(x); \
         m = 6 - m;          \
         x = x + dx );       \
       int = dx*int/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
       int result;

* now we use a known integral to test the macro
* integral of sin(x) for [0,1] = 1 - cos(1);

scalar res,diff; res = 1-cos(1); display res;

simpson(sin,0,1,2);    diff = abs(int-res); display diff;
simpson(sin,0,1,10);   diff = abs(int-res); display diff;
simpson(sin,0,1,1000); diff = abs(int-res); display diff;

* now we define a macro for the function
$macro g(y) sin(y)

simpson(g,0,1,10);    diff = abs(int-res); 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) = int;
);
display result;
gams/an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule.txt · Last modified: 2009/08/24 10:04 by support