CLHEP  2.4.7.2
C++ Class Library for High Energy Physics
ExtendedButcherTableau.icc
Go to the documentation of this file.
1 #include <ostream> // for std::endl
2 namespace Genfun {
3  ExtendedButcherTableau::ExtendedButcherTableau(const std::string &mname,
4  unsigned int xorder,
5  unsigned int xorderHat):_name(mname),
6  _order(xorder),
7  _orderHat(xorderHat){
8  }
9 
10 
11  const std::string & ExtendedButcherTableau::name() const {
12  return _name;
13  }
14 
15 
16 
17  unsigned int ExtendedButcherTableau::order() const{
18  return _order;
19  }
20 
21  unsigned int ExtendedButcherTableau::orderHat() const{
22  return _orderHat;
23  }
24 
25 
26 
27  unsigned int ExtendedButcherTableau::nSteps() const{
28  return (unsigned int)_A.size();
29  }
30 
31 
32 // don't generate warnings about intentional shadowing
33 #if defined __GNUC__
34  #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
35  #pragma GCC diagnostic push
36  #pragma GCC diagnostic ignored "-Wshadow"
37  #endif
38  #if __GNUC__ > 4
39  #pragma GCC diagnostic push
40  #pragma GCC diagnostic ignored "-Wshadow"
41  #endif
42 #endif
43 #if defined __INTEL_COMPILER
44  #pragma warning push
45  #pragma warning disable 1599
46 #endif
47 #ifdef __clang__
48  #pragma clang diagnostic push
49  #pragma clang diagnostic ignored "-Wshadow"
50 #endif
51 
52  double & ExtendedButcherTableau::A(unsigned int i, unsigned int j) {
53 
54  if (i>=(unsigned int)_A.size()) {
55  unsigned int newSize=i+1; // (shadowed)
56  for (unsigned long ii=0;ii<_A.size();ii++) {
57  _A[ii].resize(newSize,0.0);
58  }
59  for (unsigned int ii=(unsigned int)_A.size();ii<newSize;ii++) {
60  _A.push_back(std::vector<double>(newSize,0));
61  }
62 
63  if (j>=(unsigned int)_A[i].size()) {
64  unsigned int newSize=j+1; // (shadow)
65  for (unsigned long ii=0;ii<_A.size();ii++) {
66  _A[ii].resize(newSize,0.0);
67  }
68  }
69  }
70  return _A[i][j];
71  }
72 
73 #if defined __GNUC__
74  #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
75  #pragma GCC diagnostic pop
76  #endif
77  #if __GNUC__ > 4
78  #pragma GCC diagnostic pop
79  #endif
80 #endif
81 #if defined __INTEL_COMPILER
82  #pragma warning pop
83 #endif
84 #ifdef __clang__
85  #pragma clang diagnostic pop
86 #endif
87 
88  double & ExtendedButcherTableau::b(unsigned int i){
89  if (i>=(unsigned int)_b.size()) _b.resize(i+1);
90  return _b[i];
91  }
92 
93  double & ExtendedButcherTableau::bHat(unsigned int i){
94  if (i>=(unsigned int)_bHat.size()) _bHat.resize(i+1);
95  return _bHat[i];
96  }
97 
98  double & ExtendedButcherTableau::c(unsigned int i){
99  if (i>=(unsigned int)_c.size()) _c.resize(i+1);
100  return _c[i];
101  }
102 
103  const double & ExtendedButcherTableau::A(unsigned int i, unsigned int j) const{
104  return _A[i][j];
105  }
106 
107  const double & ExtendedButcherTableau::b(unsigned int i) const{
108  return _b[i];
109  }
110 
111  const double & ExtendedButcherTableau::bHat(unsigned int i) const{
112  return _bHat[i];
113  }
114 
115  const double & ExtendedButcherTableau::c(unsigned int i) const{
116  return _c[i];
117  }
118 }
119 
120 std::ostream & operator << (std::ostream & o, const Genfun::ExtendedButcherTableau & b) {
121  o << "Name " << b.name() << " of order " << b.order() << std::endl;
122  o << "A" << std::endl;
123  for (unsigned int i=0;i<b.nSteps();i++) {
124  for (unsigned int j=0;j<b.nSteps();j++) {
125  o << b.A(i,j) << " ";
126  }
127  o << std::endl;
128  }
129  o<< std::endl;
130  o << "c" << std::endl;
131  for (unsigned int j=0;j<b.nSteps();j++) {
132  o << b.c(j) << std::endl;
133  }
134  o<< std::endl;
135  o << "b" << std::endl;
136  for (unsigned int j=0;j<b.nSteps();j++) {
137  o << b.b(j) << " ";
138  }
139  o<< std::endl;
140  o << "bHat" << std::endl;
141  for (unsigned int j=0;j<b.nSteps();j++) {
142  o << b.bHat(j) << " ";
143  }
144  o << std::endl;
145  return o;
146 }
147 
148 namespace Genfun {
149 
150 
151 
152  HeunEulerXtTableau::HeunEulerXtTableau():
153  ExtendedButcherTableau("Heun-Euler Embedded Method", 2,1)
154  {
155  A(0,0) ; A(0,1) ;
156  A(1,0)=1 ; A(1,1) ;
157 
158  c(0)=0;
159  c(1)=1;
160 
161  b(0)=1/2.0;
162  b(1)=1/2.0;
163 
164  bHat(0)=1;
165  bHat(1)=0;
166  }
167 
168 
169 
170  BogackiShampineXtTableau::BogackiShampineXtTableau():
171  ExtendedButcherTableau("Bogacki-Shampine Embedded Method", 3,2)
172  {
173  A(0,0); A(0,1); A(0,2); A(0,3);
174  A(1,0)=1/2.0; A(1,1); A(1,2); A(1,3);
175  A(2,0)=0; A(2,1)=3/4.0; A(2,2); A(2,3);
176  A(3,0)=2/9.0; A(3,1)=1/3.0; A(3,2)=4/9.0; A(3,3);
177 
178  c(0)=0;
179  c(1)=1/2.;
180  c(2)=3/4.;
181  c(3)=1.0;
182 
183  b(0)=2/9.0;
184  b(1)=1/3.0;
185  b(2)=4/9.0;
186  b(3)=0;
187 
188  bHat(0)=7/24.0;
189  bHat(1)=1/4.0;
190  bHat(2)=1/3.0;
191  bHat(3)=1/8.0;
192  }
193  FehlbergRK45F2XtTableau::FehlbergRK45F2XtTableau():
194  ExtendedButcherTableau("FehlbergRK4(5) method formula 2", 4, 5)
195  {
196  A(0,0) ; A(0,1) ; A(0,2) ; A(0,3) ; A(0,4) ; A(0,5);
197  A(1,0)=1/4. ; A(1,1) ; A(1,2) ; A(1,3) ; A(1,4) ; A(1,5);
198  A(2,0)=3/32. ; A(2,1)=9/32. ; A(2,2) ; A(2,3) ; A(2,4) ; A(2,5);
199  A(3,0)=1932/2197. ; A(3,1)=-7020/2197. ; A(3,2)=7296/2197. ; A(3,3) ; A(3,4) ; A(3,5);
200  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);
201  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);
202 
203  c(0)=0;
204  c(1)=1/4.;
205  c(2)=3/8.;
206  c(3)=12/13.;
207  c(4)=1;
208  c(5)=1/2.;
209 
210  b(0)=25/216.;
211  b(1)=0;
212  b(2)=1408/2565.;
213  b(3)=2197/4104.;
214  b(4)=-1/5.;
215  b(5)=0;
216 
217  bHat(0)=16/135.;
218  bHat(1)=0;
219  bHat(2)=6656/12825.;
220  bHat(3)=28561/56430.;
221  bHat(4)=-9/50.;
222  bHat(5)=2/55.;
223 
224  }
225 
226  CashKarpXtTableau::CashKarpXtTableau():
227  ExtendedButcherTableau("FehlbergRK4(5) method formula 2", 4, 5) {
228 
229  A(0,0) ; A(0,1) ; A(0,2) ; A(0,3) ; A(0,4) ; A(0,5);
230  A(1,0) = 1/5. ; A(1,1) ; A(1,2) ; A(1,3) ; A(1,4) ; A(1,5);
231  A(2,0) = 3/40. ; A(2,1)=9/40. ; A(2,2) ; A(2,3) ; A(2,4) ; A(2,5);
232  A(3,0) = 3/10. ; A(3,1)=-9/10. ; A(3,2)=6/5. ; A(3,3) ; A(3,4) ; A(3,5);
233  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);
234  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);
235 
236  c(0)=0;
237  c(1)=1/5.0;
238  c(2)=3/10.;
239  c(3)=3/5.;
240  c(4)=1;
241  c(5)=7/8.;
242 
243  b(0)=37/378.;
244  b(1)=0;
245  b(2)=250/621.;
246  b(3)=125/594.;
247  b(4)=0;
248  b(5)=512/1771.;
249 
250  bHat(0)=2825/27648.;
251  bHat(1)=0;
252  bHat(2)=18575/48384.;
253  bHat(3)=13525/55296.;
254  bHat(4)=277/14336.;
255  bHat(5)=1/4.;
256 
257  }
258 
259 }
260 
operator<<
std::ostream & operator<<(std::ostream &o, const Genfun::ExtendedButcherTableau &b)
Definition: ExtendedButcherTableau.icc:120
Genfun
Definition: Airy.icc:9
Generated by   doxygen 1.8.14