User Tools

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.

Link to this comparison view

Both sides previous revision Previous revision
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
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;
 </​code>​ </​code>​
 +
 +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]].
IMPRESSUM / LEGAL NOTICEPRIVACY POLICY gams/an_example_how_to_use_the_gams_macro_feature_to_implement_the_simpson_s_rule.txt · Last modified: 2020/05/28 17:44 by Michael Bussieck