CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
ExtendedButcherTableau.icc
Go to the documentation of this file.
1#include <ostream> // for std::endl
2namespace Genfun {
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
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
120std::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
148namespace Genfun {
149
150
151
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
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 }
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
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
std::ostream & operator<<(std::ostream &o, const Genfun::ExtendedButcherTableau &b)
ExtendedButcherTableau(const std::string &name, unsigned int order, unsigned int orderHat)
double & A(unsigned int i, unsigned int j)
Definition Abs.hh:14