# GAMS Support Wiki

### Site Tools

gams:an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule

# Differences

This shows you the differences between two versions of the page.

 gams:an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule [2009/08/24 12:04]support gams:an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule [2020/05/28 17:44] (current)Michael Bussieck Both sides previous revision Previous revision 2020/05/28 17:44 Michael Bussieck 2009/08/24 12:04 support 2009/08/24 12:02 support created 2020/05/28 17:44 Michael Bussieck 2009/08/24 12:04 support 2009/08/24 12:02 support created Line 1: Line 1: - ====== ​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 [[http://​en.wikipedia.org/​wiki/​Simpson%27s_rule|Simpson'​s rule]] to calculate the value of a definite integral, One can even you the GAMS [[https://​www.gams.com/​latest/​docs/​UG_DollarControlOptions.html#​INDEX_macro|macro facility]] to make this look decent: <​code>​ <​code>​ - * Composit ​Simpson rule - for details see Simpson'​s rule at Wikipedia (http://​en.wikipedia.org/​wiki/​Simpson%27s_rule) + * Composite ​Simpson rule - \$macro simpson(f,​a,​b,​nn) ​    ​\ + \$macro simpson(f,​a,​b,​nn) ​                                           \ - n  = max(2,ceil(nn)); ​\ + abort\$(mod(nn,2) or nn<>​round(nn)) 'nn not an even integer'​; \ - n = n + mod(n,​2); ​    \ + dx = (b-a)/nn;                                               ​\ - dx = (b-a)/n;         ​\ + x = a + dx;                                                  \ - x = a + dx;           ​\ + interval ​= f(a) + f(b);                                      \ - int = f(a) + f(b);    \ + m = 4;                                                       ​\ - m = 4;                \ + ​for(ii = 2 to nn,                                            \ - ​for(ii = 2 to n,      \ + interval ​= interval ​+ m*f(x); ​                             \ - int = int + m*f(x); \ + m = 6 - m;                                                 ​\ - m = 6 - m;          \ + x = x + dx );                                              \ - x = x + dx );       ​\ + interval ​= dx*interval/3; - int = dx*int/3;       \ + * following will be used by the macro simpson * following will be used by the macro simpson - scalar n  number of intervals - has to be even + scalar n        '​number of intervals - has to be even' - ​a ​ interval low + ​a ​       '​interval low' - ​b ​ interval high + ​b ​       '​interval high' - dx sub interval + ​dx ​      'sub interval' - ​m ​ constant + ​m ​       '​constant' - ii loop index + ​ii ​      'loop index' - ​x ​ function argument + ​x ​       '​function argument' - int result; + interval 'result'; * now we use a known integral to test the macro * now we use a known integral to test the macro * integral of sin(x) for [0,1] = 1 - cos(1); * integral of sin(x) for [0,1] = 1 - cos(1); - scalar ​res,​diff; ​res = 1-cos(1); display ​res; + scalar ​exactResult, diff; exactResult ​= 1-cos(1); display ​exactResult; - simpson(sin,​0,​1,​2); ​   diff = abs(int-res); display diff; + simpson(sin,​0,​1,​2); ​   diff = abs(interval-exactResult); display diff; - simpson(sin,​0,​1,​10); ​  diff = abs(int-res); display diff; + simpson(sin,​0,​1,​10); ​  diff = abs(interval-exactResult); display diff; - simpson(sin,​0,​1,​1000);​ diff = abs(int-res); display diff; + simpson(sin,​0,​1,​1000);​ diff = abs(interval-exactResult); display diff; * now we define a macro for the function * now we define a macro for the function \$macro g(y) sin(y) \$macro g(y) sin(y) - simpson(g,​0,​1,​10); ​   diff = abs(int-res); display diff; + simpson(g,​0,​1,​10); ​   diff = abs(interval-exactResult); display diff; * and an even more complicated one - normal density * and an even more complicated one - normal density \$macro p(x) 1/​(sigma*sqrt(2*PI))*exp(-sqr(x-mu)/​(2*sqr(sigma))) \$macro p(x) 1/​(sigma*sqrt(2*PI))*exp(-sqr(x-mu)/​(2*sqr(sigma))) - parameters sigma,mu; + parameters sigma, mu; - set i /1*20 /, + set i / 1*20 /, - j /1*25 /; + j / 1*25 /; - parameters sx(i),​mx(i,​j),​result(i,​j);​ + parameters sx(i), mx(i,j), result(i,​j);​ sx(i) = uniform(1,​2);​ mx(i,j) = uniform(-1,​+1);​ sx(i) = uniform(1,​2);​ mx(i,j) = uniform(-1,​+1);​ Line 51: Line 53: sigma = sx(i); mu = mx(i,j); sigma = sx(i); mu = mx(i,j); ​simpson(p,​0,​10,​100);​ ​simpson(p,​0,​10,​100);​ - ​result(i,​j) = int; + ​result(i,​j) = interval; ); ); display result; display result; ​ + + This implementation has the disadvantage that it can't be used in model equations. GAMS has the ability to use [[https://​www.gams.com/​latest/​docs/​UG_ExtrinsicFunctions.html|extrinsic functions]] and [[https://​www.gams.com/​latest/​docs/​UG_ExternalEquations.html|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 [[https://​www.gams.com/​latest/​docs/​UG_ExtrinsicFunctions.html#​UG_ExtrinsicFunctions_ExternalEquations|example of an integral]].