CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
ButcherTableau.icc
Go to the documentation of this file.
1#include <ostream> // for std::endl
2namespace Genfun {
3 ButcherTableau::ButcherTableau(const std::string &xname, unsigned int xorder):_name(xname),_order(xorder){
4 }
5
6
7 const std::string & ButcherTableau::name() const {
8 return _name;
9 }
10
11
12
13 unsigned int ButcherTableau::order() const{
14 return _order;
15 }
16
17
18
19 unsigned int ButcherTableau::nSteps() const{
20 return (unsigned int)_A.size();
21 }
22
23// don't generate warnings about intentional shadowing
24#if defined __GNUC__
25 #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
26 #pragma GCC diagnostic push
27 #pragma GCC diagnostic ignored "-Wshadow"
28 #endif
29 #if __GNUC__ > 4
30 #pragma GCC diagnostic push
31 #pragma GCC diagnostic ignored "-Wshadow"
32 #endif
33#endif
34#if defined __INTEL_COMPILER
35 #pragma warning push
36 #pragma warning disable 1599
37#endif
38#ifdef __clang__
39 #pragma clang diagnostic push
40 #pragma clang diagnostic ignored "-Wshadow"
41#endif
42
43 double & ButcherTableau::A(unsigned int i, unsigned int j) {
44
45 if (i>=(unsigned int)_A.size()) {
46 unsigned int newSize=i+1; // (shadowed)
47 for (unsigned long ii=0;ii<_A.size();ii++) {
48 _A[ii].resize(newSize,0.0);
49 }
50 for (unsigned int ii=(unsigned int)_A.size();ii<newSize;ii++) {
51 _A.push_back(std::vector<double>(newSize,0));
52 }
53
54 if (j>=(unsigned int)_A[i].size()) {
55 unsigned int newSize=j+1; // (shadow)
56 for (unsigned long ii=0;ii<_A.size();ii++) {
57 _A[ii].resize(newSize,0.0);
58 }
59 }
60 }
61 return _A[i][j];
62 }
63
64#if defined __GNUC__
65 #if __GNUC__ > 3 && __GNUC_MINOR__ > 6
66 #pragma GCC diagnostic pop
67 #endif
68 #if __GNUC__ > 4
69 #pragma GCC diagnostic pop
70 #endif
71#endif
72#if defined __INTEL_COMPILER
73 #pragma warning pop
74#endif
75#ifdef __clang__
76 #pragma clang diagnostic pop
77#endif
78
79 double & ButcherTableau::b(unsigned int i){
80 if (i>=(unsigned int)_b.size()) _b.resize(i+1);
81 return _b[i];
82 }
83
84 double & ButcherTableau::c(unsigned int i){
85 if (i>=(unsigned int)_c.size()) _c.resize(i+1);
86 return _c[i];
87 }
88
89 const double & ButcherTableau::A(unsigned int i, unsigned int j) const{
90 return _A[i][j];
91 }
92
93 const double & ButcherTableau::b(unsigned int i) const{
94 return _b[i];
95 }
96
97 const double & ButcherTableau::c(unsigned int i) const{
98 return _c[i];
99 }
100}
101
102std::ostream & operator << (std::ostream & o, const Genfun::ButcherTableau & b) {
103 o << "Name " << b.name() << " of order " << b.order() << std::endl;
104 o << "A" << std::endl;
105 for (unsigned int i=0;i<b.nSteps();i++) {
106 for (unsigned int j=0;j<b.nSteps();j++) {
107 o << b.A(i,j) << " ";
108 }
109 o << std::endl;
110 }
111 o<< std::endl;
112 o << "c" << std::endl;
113 for (unsigned int j=0;j<b.nSteps();j++) {
114 o << b.c(j) << std::endl;
115 }
116 o<< std::endl;
117 o << "b" << std::endl;
118 for (unsigned int j=0;j<b.nSteps();j++) {
119 o << b.b(j) << " ";
120 }
121 o << std::endl;
122 return o;
123}
124
125namespace Genfun {
127 ButcherTableau("Euler Method", 1)
128 {
129 A(0,0)=0;
130 b(0)=1;
131 c(0)=1;
132 }
133
135 ButcherTableau("Midpoint Method", 2)
136 {
137 A(1,0)=1/2.0;
138 c(0)=0;
139 c(1)=1/2.0;
140 b(0)=0;
141 b(1)=1;
142
143 }
144
146 ButcherTableau("Trapezoid Method", 2)
147 {
148 A(1,0)=1;
149 c(0)=0;
150 c(1)=1;
151 b(0)=1/2.0;
152 b(1)=1/2.0;
153
154 }
155
157 ButcherTableau("RK31 Method", 3)
158 {
159 A(0,0) ; A(0,1) ; A(0,2);
160 A(1,0)=2/3.0; A(1,1) ; A(1,2);
161 A(2,0)=1/3.0; A(2,1)=1/3.0; A(2,2);
162
163 c(0)=0;
164 c(1)=2/3.0;
165 c(2)=2/3.0;
166 b(0)=1/4.0;
167 b(1)=0;
168 b(2)=3/4.0;
169 }
170
171
173 ButcherTableau("RK32 Method", 3)
174 {
175 A(0,0) ; A(0,1) ; A(0,2);
176 A(1,0)=1/2.0; A(1,1) ; A(1,2);
177 A(2,0)=-1 ; A(2,1)= 2 ; A(2,2);
178
179 c(0)=0;
180 c(1)=1/2.0;
181 c(2)=1;
182 b(0)=1/6.0;
183 b(1)=2/3.0;
184 b(2)=1/6.0;
185
186 }
187
189 ButcherTableau("Classical Runge Kutta Method", 4)
190 {
191 A(0,0) ; A(0,1) ; A(0,2) ; A(0,3);
192 A(1,0)=1/2.0; A(1,1) ; A(1,2) ; A(1,3);
193 A(2,0)=0 ; A(2,1)=1/2.0 ; A(2,2) ; A(2,3);
194 A(3,0)=0 ; A(3,1)=0 ; A(3,2)=1 ; A(3,3);
195
196 c(0)=0;
197 c(1)=1/2.0;
198 c(2)=1/2.0;
199 c(3)=1;
200 b(0)=1/6.0;
201 b(1)=1/3.0;
202 b(2)=1/3.0;
203 b(3)=1/6.0;
204 }
205
207 ButcherTableau("Three-Eighths Rule Method", 4)
208 {
209 A(0,0) ; A(0,1) ; A(0,2) ; A(0,3);
210 A(1,0)=1/3.0 ; A(1,1) ; A(1,2) ; A(1,3);
211 A(2,0)=-1/3.0 ; A(2,1)=1 ; A(2,2) ; A(2,3);
212 A(3,0)=1 ; A(3,1)=-1 ; A(3,2)=1 ; A(3,3);
213
214 c(0)=0;
215 c(1)=1/3.0;
216 c(2)=2/3.0;
217 c(3)=1;
218 b(0)=1/8.0;
219 b(1)=3/8.0;
220 b(2)=3/8.0;
221 b(3)=1/8.0;
222 }
223}
224
std::ostream & operator<<(std::ostream &o, const Genfun::ButcherTableau &b)
ButcherTableau(const std::string &name, unsigned int order)
unsigned int nSteps() const
const std::string & name() const
double & A(unsigned int i, unsigned int j)
double & c(unsigned int i)
double & b(unsigned int i)
unsigned int order() const
Definition Abs.hh:14