Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:54:34

0001 #include <ostream> // for std::endl
0002 namespace Genfun {
0003   ExtendedButcherTableau::ExtendedButcherTableau(const std::string &mname, 
0004                          unsigned int xorder,
0005                          unsigned int xorderHat):_name(mname),
0006                                     _order(xorder),
0007                                     _orderHat(xorderHat){
0008   }
0009   
0010   
0011   const std::string & ExtendedButcherTableau::name() const {
0012     return _name;
0013   }
0014   
0015   
0016   
0017   unsigned int ExtendedButcherTableau::order() const{
0018     return _order;
0019   }
0020 
0021   unsigned int ExtendedButcherTableau::orderHat() const{
0022     return _orderHat;
0023   }
0024   
0025   
0026   
0027   unsigned int ExtendedButcherTableau::nSteps() const{
0028     return (unsigned int)_A.size();
0029   }
0030   
0031 
0032 // don't generate warnings about intentional shadowing
0033 #if defined __GNUC__ 
0034   #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
0035     #pragma GCC diagnostic push
0036     #pragma GCC diagnostic ignored "-Wshadow"
0037   #endif
0038   #if __GNUC__ > 4
0039     #pragma GCC diagnostic push
0040     #pragma GCC diagnostic ignored "-Wshadow"
0041   #endif
0042 #endif
0043 #if defined  __INTEL_COMPILER
0044   #pragma warning push
0045   #pragma warning disable 1599
0046 #endif
0047 #ifdef __clang__
0048   #pragma clang diagnostic push
0049   #pragma clang diagnostic ignored "-Wshadow"
0050 #endif
0051   
0052   double & ExtendedButcherTableau::A(unsigned int i, unsigned int j) {
0053     
0054     if (i>=(unsigned int)_A.size()) {
0055       unsigned int newSize=i+1; // (shadowed)
0056       for (unsigned long ii=0;ii<_A.size();ii++) {
0057     _A[ii].resize(newSize,0.0);
0058       }
0059       for (unsigned int ii=(unsigned int)_A.size();ii<newSize;ii++) {
0060     _A.push_back(std::vector<double>(newSize,0));
0061       }
0062       
0063       if (j>=(unsigned int)_A[i].size()) {
0064     unsigned int newSize=j+1;  // (shadow)
0065     for (unsigned long ii=0;ii<_A.size();ii++) {
0066       _A[ii].resize(newSize,0.0);
0067     }
0068       }
0069     }
0070     return _A[i][j];
0071   }
0072 
0073 #if defined __GNUC__ 
0074   #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
0075     #pragma GCC diagnostic pop
0076   #endif
0077   #if __GNUC__ > 4
0078     #pragma GCC diagnostic pop
0079   #endif
0080 #endif
0081 #if defined  __INTEL_COMPILER
0082   #pragma warning pop
0083 #endif
0084 #ifdef __clang__
0085   #pragma clang diagnostic pop
0086 #endif
0087 
0088   double & ExtendedButcherTableau::b(unsigned int i){
0089     if (i>=(unsigned int)_b.size()) _b.resize(i+1);
0090     return _b[i];
0091   }
0092 
0093   double & ExtendedButcherTableau::bHat(unsigned int i){
0094     if (i>=(unsigned int)_bHat.size()) _bHat.resize(i+1);
0095     return _bHat[i];
0096   }
0097 
0098   double & ExtendedButcherTableau::c(unsigned int i){
0099     if (i>=(unsigned int)_c.size()) _c.resize(i+1);
0100     return _c[i];
0101   }
0102   
0103   const double & ExtendedButcherTableau::A(unsigned int i, unsigned int j) const{
0104     return _A[i][j];
0105   }
0106   
0107   const double & ExtendedButcherTableau::b(unsigned int i) const{
0108     return _b[i];
0109   }
0110   
0111   const double & ExtendedButcherTableau::bHat(unsigned int i) const{
0112     return _bHat[i];
0113   }
0114   
0115   const double & ExtendedButcherTableau::c(unsigned int i) const{
0116     return _c[i];
0117   }
0118 }
0119 
0120 std::ostream & operator << (std::ostream & o, const Genfun::ExtendedButcherTableau & b) {
0121   o << "Name " << b.name() << " of order " << b.order() << std::endl;
0122   o << "A" << std::endl;
0123   for (unsigned int i=0;i<b.nSteps();i++) {
0124     for (unsigned int j=0;j<b.nSteps();j++) {
0125       o << b.A(i,j) << " ";
0126     }
0127     o << std::endl;
0128   }
0129   o<< std::endl;
0130   o << "c" << std::endl;
0131   for (unsigned int j=0;j<b.nSteps();j++) {
0132     o << b.c(j) << std::endl;
0133   }
0134   o<< std::endl;
0135   o << "b" << std::endl;
0136   for (unsigned int j=0;j<b.nSteps();j++) {
0137     o << b.b(j) << " ";
0138   }
0139   o<< std::endl;
0140   o << "bHat" << std::endl;
0141   for (unsigned int j=0;j<b.nSteps();j++) {
0142     o << b.bHat(j) << " ";
0143   }
0144   o << std::endl;
0145   return o;
0146 }
0147 
0148 namespace Genfun {
0149 
0150 
0151 
0152   HeunEulerXtTableau::HeunEulerXtTableau():
0153     ExtendedButcherTableau("Heun-Euler Embedded Method", 2,1)
0154   {
0155     A(0,0)      ;    A(0,1)      ;
0156     A(1,0)=1    ;    A(1,1)      ;
0157 
0158     c(0)=0;
0159     c(1)=1;
0160 
0161     b(0)=1/2.0;
0162     b(1)=1/2.0;
0163 
0164     bHat(0)=1;
0165     bHat(1)=0;
0166   }
0167 
0168 
0169 
0170   BogackiShampineXtTableau::BogackiShampineXtTableau():
0171     ExtendedButcherTableau("Bogacki-Shampine Embedded Method", 3,2)
0172   {
0173     A(0,0);         A(0,1);        A(0,2);       A(0,3);
0174     A(1,0)=1/2.0;   A(1,1);        A(1,2);       A(1,3);
0175     A(2,0)=0;       A(2,1)=3/4.0;  A(2,2);       A(2,3);
0176     A(3,0)=2/9.0;   A(3,1)=1/3.0;  A(3,2)=4/9.0; A(3,3);
0177 
0178     c(0)=0;
0179     c(1)=1/2.;
0180     c(2)=3/4.;
0181     c(3)=1.0;
0182 
0183     b(0)=2/9.0;
0184     b(1)=1/3.0;
0185     b(2)=4/9.0;
0186     b(3)=0;
0187 
0188     bHat(0)=7/24.0;
0189     bHat(1)=1/4.0;
0190     bHat(2)=1/3.0;
0191     bHat(3)=1/8.0;
0192   }
0193   FehlbergRK45F2XtTableau::FehlbergRK45F2XtTableau():
0194   ExtendedButcherTableau("FehlbergRK4(5) method formula 2", 4, 5)
0195   {
0196     A(0,0)               ;   A(0,1)                     ;        A(0,2)                 ;       A(0,3)           ;     A(0,4)        ;   A(0,5);
0197     A(1,0)=1/4.          ;   A(1,1)                     ;        A(1,2)                 ;       A(1,3)           ;     A(1,4)        ;   A(1,5);
0198     A(2,0)=3/32.         ;   A(2,1)=9/32.               ;        A(2,2)                 ;       A(2,3)           ;     A(2,4)        ;   A(2,5);
0199     A(3,0)=1932/2197.    ;   A(3,1)=-7020/2197.         ;        A(3,2)=7296/2197.      ;       A(3,3)           ;     A(3,4)        ;   A(3,5);
0200     A(4,0)=439/216.      ;   A(4,1)=-8.                 ;        A(4,2)=3680/513.       ;       A(4,3)=-845/4104.;     A(4,4)        ;   A(4,5);
0201     A(5,0)=-8/27.        ;   A(5,1)=2.                  ;        A(5,2)=-3544/2565.     ;       A(5,3)=1859/4104.;     A(5,4)=-11/40.;   A(5,5);
0202 
0203     c(0)=0;
0204     c(1)=1/4.;
0205     c(2)=3/8.;
0206     c(3)=12/13.;
0207     c(4)=1;
0208     c(5)=1/2.;
0209 
0210     b(0)=25/216.;
0211     b(1)=0;
0212     b(2)=1408/2565.;
0213     b(3)=2197/4104.;
0214     b(4)=-1/5.;
0215     b(5)=0;
0216 
0217     bHat(0)=16/135.;
0218     bHat(1)=0;
0219     bHat(2)=6656/12825.;
0220     bHat(3)=28561/56430.;
0221     bHat(4)=-9/50.;
0222     bHat(5)=2/55.;
0223 
0224   }
0225 
0226   CashKarpXtTableau::CashKarpXtTableau():
0227   ExtendedButcherTableau("FehlbergRK4(5) method formula 2", 4, 5) {
0228     
0229     A(0,0)              ;   A(0,1)         ;        A(0,2)           ;       A(0,3)              ;     A(0,4)          ;  A(0,5);
0230     A(1,0) = 1/5.       ;   A(1,1)         ;        A(1,2)           ;       A(1,3)              ;     A(1,4)          ;  A(1,5);
0231     A(2,0) = 3/40.      ;   A(2,1)=9/40.   ;        A(2,2)           ;       A(2,3)              ;     A(2,4)          ;  A(2,5);
0232     A(3,0) = 3/10.      ;   A(3,1)=-9/10.  ;        A(3,2)=6/5.      ;       A(3,3)              ;     A(3,4)          ;  A(3,5);
0233     A(4,0) = -11/54.    ;   A(4,1)=5/2.    ;        A(4,2)=-70/27.   ;       A(4,3)=35/27.       ;     A(4,4)          ;  A(4,5);
0234     A(5,0) = 1631/55296.;   A(5,1)=175/512.;        A(5,2)=575/13824.;       A(5,3)=44275/110592.;     A(5,4)=253/4096.;  A(5,5);
0235 
0236     c(0)=0;
0237     c(1)=1/5.0;
0238     c(2)=3/10.;
0239     c(3)=3/5.;
0240     c(4)=1;
0241     c(5)=7/8.;
0242 
0243     b(0)=37/378.;
0244     b(1)=0;
0245     b(2)=250/621.;
0246     b(3)=125/594.;
0247     b(4)=0;
0248     b(5)=512/1771.;
0249 
0250     bHat(0)=2825/27648.;
0251     bHat(1)=0;
0252     bHat(2)=18575/48384.;
0253     bHat(3)=13525/55296.;
0254     bHat(4)=277/14336.;
0255     bHat(5)=1/4.;
0256 
0257   }
0258 
0259 }
0260