CLHEP/GenericFunctions/ButcherTableau.icc

0001 #include <ostream> // for std::endl
0002 namespace Genfun {
0003   ButcherTableau::ButcherTableau(const std::string &xname, unsigned int xorder):_name(xname),_order(xorder){
0004   }
0005
0006
0007   const std::string & ButcherTableau::name() const {
0008     return _name;
0009   }
0010
0011
0012
0013   unsigned int ButcherTableau::order() const{
0014     return _order;
0015   }
0016
0017
0018
0019   unsigned int ButcherTableau::nSteps() const{
0020     return (unsigned int)_A.size();
0021   }
0022
0023 // don't generate warnings about intentional shadowing
0024 #if defined __GNUC__
0025   #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
0026     #pragma GCC diagnostic push
0027     #pragma GCC diagnostic ignored "-Wshadow"
0028   #endif
0029   #if __GNUC__ > 4
0030     #pragma GCC diagnostic push
0031     #pragma GCC diagnostic ignored "-Wshadow"
0032   #endif
0033 #endif
0034 #if defined  __INTEL_COMPILER
0035   #pragma warning push
0036   #pragma warning disable 1599
0037 #endif
0038 #ifdef __clang__
0039   #pragma clang diagnostic push
0040   #pragma clang diagnostic ignored "-Wshadow"
0041 #endif
0042
0043   double & ButcherTableau::A(unsigned int i, unsigned int j) {
0044
0045     if (i>=(unsigned int)_A.size()) {
0046       unsigned int newSize=i+1; // (shadowed)
0047       for (unsigned long ii=0;ii<_A.size();ii++) {
0048     _A[ii].resize(newSize,0.0);
0049       }
0050       for (unsigned int ii=(unsigned int)_A.size();ii<newSize;ii++) {
0051     _A.push_back(std::vector<double>(newSize,0));
0052       }
0053
0054       if (j>=(unsigned int)_A[i].size()) {
0055     unsigned int newSize=j+1;  // (shadow)
0056     for (unsigned long ii=0;ii<_A.size();ii++) {
0057       _A[ii].resize(newSize,0.0);
0058     }
0059       }
0060     }
0061     return _A[i][j];
0062   }
0063
0064 #if defined __GNUC__
0065   #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
0066     #pragma GCC diagnostic pop
0067   #endif
0068   #if __GNUC__ > 4
0069     #pragma GCC diagnostic pop
0070   #endif
0071 #endif
0072 #if defined  __INTEL_COMPILER
0073   #pragma warning pop
0074 #endif
0075 #ifdef __clang__
0076   #pragma clang diagnostic pop
0077 #endif
0078
0079   double & ButcherTableau::b(unsigned int i){
0080     if (i>=(unsigned int)_b.size()) _b.resize(i+1);
0081     return _b[i];
0082   }
0083
0084   double & ButcherTableau::c(unsigned int i){
0085     if (i>=(unsigned int)_c.size()) _c.resize(i+1);
0086     return _c[i];
0087   }
0088
0089   const double & ButcherTableau::A(unsigned int i, unsigned int j) const{
0090     return _A[i][j];
0091   }
0092
0093   const double & ButcherTableau::b(unsigned int i) const{
0094     return _b[i];
0095   }
0096
0097   const double & ButcherTableau::c(unsigned int i) const{
0098     return _c[i];
0099   }
0100 }
0101
0102 std::ostream & operator << (std::ostream & o, const Genfun::ButcherTableau & b) {
0103   o << "Name " << b.name() << " of order " << b.order() << std::endl;
0104   o << "A" << std::endl;
0105   for (unsigned int i=0;i<b.nSteps();i++) {
0106     for (unsigned int j=0;j<b.nSteps();j++) {
0107       o << b.A(i,j) << " ";
0108     }
0109     o << std::endl;
0110   }
0111   o<< std::endl;
0112   o << "c" << std::endl;
0113   for (unsigned int j=0;j<b.nSteps();j++) {
0114     o << b.c(j) << std::endl;
0115   }
0116   o<< std::endl;
0117   o << "b" << std::endl;
0118   for (unsigned int j=0;j<b.nSteps();j++) {
0119     o << b.b(j) << " ";
0120   }
0121   o << std::endl;
0122   return o;
0123 }
0124
0125 namespace Genfun {
0126   EulerTableau::EulerTableau():
0127     ButcherTableau("Euler Method", 1)
0128   {
0129     A(0,0)=0;
0130     b(0)=1;
0131     c(0)=1;
0132   }
0133
0134   MidpointTableau::MidpointTableau():
0135     ButcherTableau("Midpoint Method", 2)
0136   {
0137     A(1,0)=1/2.0;
0138     c(0)=0;
0139     c(1)=1/2.0;
0140     b(0)=0;
0141     b(1)=1;
0142
0143   }
0144
0145   TrapezoidTableau::TrapezoidTableau():
0146     ButcherTableau("Trapezoid Method", 2)
0147   {
0148     A(1,0)=1;
0149     c(0)=0;
0150     c(1)=1;
0151     b(0)=1/2.0;
0152     b(1)=1/2.0;
0153
0154   }
0155
0156   RK31Tableau::RK31Tableau():
0157     ButcherTableau("RK31 Method", 3)
0158   {
0159     A(0,0)      ;    A(0,1)      ;   A(0,2);
0160     A(1,0)=2/3.0;    A(1,1)      ;   A(1,2);
0161     A(2,0)=1/3.0;    A(2,1)=1/3.0;   A(2,2);
0162
0163     c(0)=0;
0164     c(1)=2/3.0;
0165     c(2)=2/3.0;
0166     b(0)=1/4.0;
0167     b(1)=0;
0168     b(2)=3/4.0;
0169   }
0170
0171
0172   RK32Tableau::RK32Tableau():
0173     ButcherTableau("RK32 Method", 3)
0174   {
0175     A(0,0)      ;    A(0,1)      ;   A(0,2);
0176     A(1,0)=1/2.0;    A(1,1)      ;   A(1,2);
0177     A(2,0)=-1   ;    A(2,1)= 2   ;   A(2,2);
0178
0179     c(0)=0;
0180     c(1)=1/2.0;
0181     c(2)=1;
0182     b(0)=1/6.0;
0183     b(1)=2/3.0;
0184     b(2)=1/6.0;
0185
0186   }
0187
0188   ClassicalRungeKuttaTableau::ClassicalRungeKuttaTableau():
0189     ButcherTableau("Classical Runge Kutta Method", 4)
0190   {
0191     A(0,0)      ;    A(0,1)       ;   A(0,2)    ; A(0,3);
0192     A(1,0)=1/2.0;    A(1,1)       ;   A(1,2)    ; A(1,3);
0193     A(2,0)=0    ;    A(2,1)=1/2.0 ;   A(2,2)    ; A(2,3);
0194     A(3,0)=0    ;    A(3,1)=0     ;   A(3,2)=1  ; A(3,3);
0195
0196     c(0)=0;
0197     c(1)=1/2.0;
0198     c(2)=1/2.0;
0199     c(3)=1;
0200     b(0)=1/6.0;
0201     b(1)=1/3.0;
0202     b(2)=1/3.0;
0203     b(3)=1/6.0;
0204   }
0205
0206   ThreeEighthsRuleTableau::ThreeEighthsRuleTableau():
0207     ButcherTableau("Three-Eighths Rule Method", 4)
0208   {
0209     A(0,0)        ;    A(0,1)       ;   A(0,2)    ; A(0,3);
0210     A(1,0)=1/3.0  ;    A(1,1)       ;   A(1,2)    ; A(1,3);
0211     A(2,0)=-1/3.0 ;    A(2,1)=1     ;   A(2,2)    ; A(2,3);
0212     A(3,0)=1      ;    A(3,1)=-1    ;   A(3,2)=1  ; A(3,3);
0213
0214     c(0)=0;
0215     c(1)=1/3.0;
0216     c(2)=2/3.0;
0217     c(3)=1;
0218     b(0)=1/8.0;
0219     b(1)=3/8.0;
0220     b(2)=3/8.0;
0221     b(3)=1/8.0;
0222   }
0223 }
0224