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
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;
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;
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