HermesCommon  3.0
table.cpp
1 // This file is part of HermesCommon
2 //
3 // Hermes is free software; you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation; either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes; if not, see <http://www.gnu.prg/licenses/>.
18 #include "table.h"
19 #include "matrix.h"
20 #include "dense_matrix_operations.h"
21 #include "util/memory_handling.h"
22 
23 using namespace Hermes::Algebra::DenseMatrixOperations;
24 
25 namespace Hermes
26 {
27  Table::Table()
28  {
29  this->size = 0;
30  this->A = nullptr;
31  }
32 
33  Table::Table(unsigned int size)
34  {
35  // Size.
36  this->size = size;
37  // A array.
38  this->A = new_matrix<double>(size, size);
39  for (unsigned int i = 0; i < size; i++)
40  {
41  for (unsigned int j = 0; j < size; j++) this->A[i][j] = 0;
42  }
43  }
44 
45  void Table::alloc(unsigned int size)
46  {
47  // Size.
48  this->size = size;
49  // A array.
50  this->A = new_matrix<double>(size, size);
51  for (unsigned int i = 0; i < size; i++)
52  {
53  for (unsigned int j = 0; j < size; j++) this->A[i][j] = 0;
54  }
55  }
56 
57  unsigned int Table::get_size()
58  {
59  return this->size;
60  }
61 
62  double Table::get_A(unsigned int i, unsigned int j)
63  {
64  if ((i >= this->size) || (j >= this->size))
65  throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
66  return this->A[i][j];
67  }
68 
69  void Table::set_A(unsigned int i, unsigned int j, double val)
70  {
71  if ((i >= this->size) || (j >= this->size))
72  throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
73  this->A[i][j] = val;
74  }
75 
76  ButcherTable::ButcherTable() : Table()
77  {
78  this->B = nullptr;
79  this->B2 = nullptr;
80  this->C = nullptr;
81  }
82 
83  ButcherTable::ButcherTable(unsigned int size) : Table(size)
84  {
85  // B array.
86  this->B = malloc_with_check<ButcherTable, double>(size, this);
87  for (unsigned int j = 0; j < size; j++) this->B[j] = 0;
88  // B2 array.
89  this->B2 = malloc_with_check<ButcherTable, double>(size, this);
90  for (unsigned int j = 0; j < size; j++) this->B2[j] = 0;
91  // C array.
92  this->C = malloc_with_check<ButcherTable, double>(size, this);
93  for (unsigned int j = 0; j < size; j++) this->C[j] = 0;
94  }
95 
96  ButcherTable::ButcherTable(ButcherTableType butcher_table)
97  {
98  switch (butcher_table)
99  {
100  /* EXPLICIT METHODS */
101 
102  // Explicit Euler.
103  case Explicit_RK_1:
104  this->alloc(1);
105  this->set_B(0, 1.);
106  break;
107 
108  // Explicit RK-2.
109  case Explicit_RK_2:
110  this->alloc(2);
111  this->set_A(1, 0, 2. / 3.);
112  this->set_A(1, 1, 0.);
113  this->set_B(0, 1. / 4.);
114  this->set_B(1, 3. / 4.);
115  this->set_C(0, 2. / 3.);
116  break;
117 
118  // Explicit RK-3.
119  case Explicit_RK_3:
120  this->alloc(3);
121  this->set_A(1, 0, 1. / 2.);
122  this->set_A(2, 0, -1.);
123  this->set_A(2, 1, 2.);
124  this->set_B(0, 1. / 6.);
125  this->set_B(1, 2. / 3.);
126  this->set_B(2, 1. / 6.);
127  this->set_C(1, 1. / 2.);
128  this->set_C(2, 1.);
129  break;
130 
131  // Explicit RK-4.
132  case Explicit_RK_4:
133  this->alloc(4);
134  this->set_A(1, 0, 1. / 2.);
135  this->set_A(2, 1, 1. / 2.);
136  this->set_A(3, 2, 1.);
137  this->set_B(0, 1. / 6.);
138  this->set_B(1, 1. / 3.);
139  this->set_B(2, 1. / 3.);
140  this->set_B(3, 1. / 6.);
141  this->set_C(1, 1. / 2.);
142  this->set_C(2, 1. / 2.);
143  this->set_C(3, 1.);
144  break;
145 
146  /* IMPLICIT METHODS */
147 
148  // Implicit Euler.
149  case Implicit_RK_1:
150  this->alloc(1);
151  this->set_A(0, 0, 1.);
152  this->set_B(0, 1.);
153  this->set_C(0, 1.);
154  break;
155 
156  // Implicit Crank Nicolson.
157  case Implicit_Crank_Nicolson_2_2:
158  this->alloc(2);
159  this->set_A(0, 0, 1. / 2.);
160  this->set_A(0, 1, 1. / 2.);
161  this->set_B(0, 1. / 2.);
162  this->set_B(1, 1. / 2.);
163  this->set_C(0, 1.);
164  break;
165 
166  // Implicit SIRK-22 (second-order).
167  case Implicit_SIRK_2_2:
168  this->alloc(2);
169  this->set_A(0, 0, (5. - 3.*sqrt(2.)) / 4.);
170  this->set_A(0, 1, (7. - 5.*sqrt(2.)) / 4.);
171  this->set_A(1, 0, (1. + 1.*sqrt(2.)) / 4.);
172  this->set_A(1, 1, (3. - 1.*sqrt(2.)) / 4.);
173  this->set_B(0, (1. + 1.*sqrt(2.)) / 4.);
174  this->set_B(1, (3. - 1.*sqrt(2.)) / 4.);
175  this->set_C(0, 3. - 2.*sqrt(2.));
176  this->set_C(1, 1.);
177  break;
178 
179  // Implicit ESIRK-22 (second-order).
180  case Implicit_ESIRK_2_2:
181  this->alloc(2);
182  this->set_A(0, 0, (9. - 6.*sqrt(2.)) / 4.);
183  this->set_A(0, 1, (-3. + 2.*sqrt(2.)) / 4.);
184  this->set_A(1, 0, (11. - 6.*sqrt(2.)) / 4.);
185  this->set_A(1, 1, (-1. + 2.*sqrt(2.)) / 4.);
186  this->set_B(0, 2. - sqrt(2.0));
187  this->set_B(1, -1. + sqrt(2.0));
188  this->set_C(1, 1.);
189  break;
190 
191  // Implicit SDIRK-2-2 (second-order).
192  case Implicit_SDIRK_2_2:
193  this->alloc(2);
194  this->set_A(0, 0, 1. - 1. / sqrt(2.));
195  this->set_A(1, 0, 1. / sqrt(2.));
196  this->set_A(1, 1, 1. - 1. / sqrt(2.));
197  this->set_B(0, 1. / sqrt(2.));
198  this->set_B(1, 1. - 1. / sqrt(2.));
199  this->set_C(0, 1. - 1. / sqrt(2.));
200  this->set_C(1, 1.);
201  break;
202 
203  // Implicit Lobatto IIIA (second-order).
204  case Implicit_Lobatto_IIIA_2_2:
205  this->alloc(2);
206  this->set_A(1, 0, 1. / 2.);
207  this->set_A(1, 1, 1. / 2.);
208  this->set_B(0, 1. / 2.);
209  this->set_B(1, 1. / 2.);
210  this->set_C(1, 1.);
211  break;
212 
213  // Implicit Lobatto IIIB (second-order).
214  case Implicit_Lobatto_IIIB_2_2:
215  this->alloc(2);
216  this->set_A(0, 0, 1. / 2.);
217  this->set_A(0, 1, 1. / 2.);
218  this->set_B(0, 1. / 2.);
219  this->set_B(1, 1. / 2.);
220  this->set_C(0, 1. / 2.);
221  this->set_C(1, 1. / 2.);
222  break;
223 
224  // Implicit Lobatto IIIC (second-order).
225  case Implicit_Lobatto_IIIC_2_2:
226  this->alloc(2);
227  this->set_A(0, 0, 1. / 2.);
228  this->set_A(0, 1, -1. / 2.);
229  this->set_A(1, 0, 1. / 2.);
230  this->set_A(1, 1, 1. / 2.);
231  this->set_B(0, 1. / 2.);
232  this->set_B(1, 1. / 2.);
233  this->set_C(1, 1.);
234  break;
235 
236  // Implicit Lobatto IIIA (fourth-order).
237  case Implicit_Lobatto_IIIA_3_4:
238  this->alloc(3);
239  this->set_A(1, 0, 5. / 24.);
240  this->set_A(2, 0, 1. / 6.);
241  this->set_A(1, 1, 1. / 3.);
242  this->set_A(2, 1, 2. / 3.);
243  this->set_A(1, 2, -1. / 24.);
244  this->set_A(2, 2, 1. / 6.);
245  this->set_B(0, 1. / 6.);
246  this->set_B(1, 2. / 3.);
247  this->set_B(2, 1. / 6.);
248  this->set_C(1, 1. / 2.);
249  this->set_C(2, 1.);
250  break;
251 
252  // Implicit Lobatto IIIB (fourth-order).
253  case Implicit_Lobatto_IIIB_3_4:
254  this->alloc(3);
255  this->set_A(0, 0, 1. / 6.);
256  this->set_A(1, 0, 1. / 6.);
257  this->set_A(2, 0, 1. / 6.);
258  this->set_A(0, 1, -1. / 6.);
259  this->set_A(1, 1, 1. / 3.);
260  this->set_A(2, 1, 5. / 6.);
261  this->set_B(0, 1. / 6.);
262  this->set_B(1, 2. / 3.);
263  this->set_B(2, 1. / 6.);
264  this->set_C(1, 1. / 2.);
265  this->set_C(2, 1.);
266  break;
267 
268  // Implicit Lobatto IIIC (fourth-order).
269  case Implicit_Lobatto_IIIC_3_4:
270  this->alloc(3);
271  this->set_A(0, 0, 1. / 6.);
272  this->set_A(1, 0, 1. / 6.);
273  this->set_A(2, 0, 1. / 6.);
274  this->set_A(0, 1, -1. / 3.);
275  this->set_A(1, 1, 5. / 12.);
276  this->set_A(2, 1, 2. / 3.);
277  this->set_A(0, 2, 1. / 6.);
278  this->set_A(1, 2, -1. / 12.);
279  this->set_A(2, 2, 1. / 6.);
280  this->set_B(0, 1. / 6.);
281  this->set_B(1, 2. / 3.);
282  this->set_B(2, 1. / 6.);
283  this->set_C(1, 1. / 2.);
284  this->set_C(2, 1.);
285  break;
286 
287  // Implicit Radau IIA (fifth-order).
288  case Implicit_Radau_IIA_3_5:
289  this->alloc(3);
290  this->set_A(0, 0, (88. - 7 * sqrt((double)6)) / 360);
291  this->set_A(1, 0, (296 + 169 * sqrt((double)6)) / 1800);
292  this->set_A(2, 0, (16. - sqrt((double)6)) / 36);
293  this->set_A(0, 1, (296 - 169 * sqrt((double)6)) / 1800);
294  this->set_A(1, 1, (88. + 7 * sqrt((double)6)) / 360);
295  this->set_A(2, 1, (16. + sqrt((double)6)) / 36);
296  this->set_A(0, 2, (-2. + 3 * sqrt((double)6)) / 225);
297  this->set_A(1, 2, (-2. - 3 * sqrt((double)6)) / 225);
298  this->set_A(2, 2, 1. / 9.);
299  this->set_B(0, (16. - sqrt((double)6)) / 36);
300  this->set_B(1, (16. + sqrt((double)6)) / 36);
301  this->set_B(2, 1. / 9.);
302  this->set_C(0, (4. - sqrt((double)6)) / 10);
303  this->set_C(1, (4. + sqrt((double)6)) / 10);
304  this->set_C(2, 1.);
305  break;
306 
307  // Implicit SDIRK-5-4 (fourth-order).
308  case Implicit_SDIRK_5_4:
309  this->alloc(5);
310  this->set_A(0, 0, 1. / 4.);
311  this->set_A(1, 0, 1. / 2.);
312  this->set_A(1, 1, 1. / 4.);
313  this->set_A(2, 0, 17. / 50.);
314  this->set_A(2, 1, -1. / 25.);
315  this->set_A(2, 2, 1. / 4.);
316  this->set_A(3, 0, 371. / 1360.);
317  this->set_A(3, 1, -137. / 2720.);
318  this->set_A(3, 2, 15. / 544.);
319  this->set_A(3, 3, 1. / 4.);
320  this->set_A(4, 0, 25. / 24.);
321  this->set_A(4, 1, -49. / 48.);
322  this->set_A(4, 2, 125. / 16.);
323  this->set_A(4, 3, -85. / 12.);
324  this->set_A(4, 4, 1. / 4.);
325  this->set_B(0, 25. / 24.);
326  this->set_B(1, -49. / 48.);
327  this->set_B(2, 125. / 16.);
328  this->set_B(3, -85. / 12.);
329  this->set_B(4, 1. / 4.);
330  this->set_C(0, 1. / 4.);
331  this->set_C(1, 3. / 4.);
332  this->set_C(2, 11. / 20.);
333  this->set_C(3, 1. / 2.);
334  this->set_C(4, 1.);
335  break;
336 
337  /* EMBEDDED EXPLICIT METHODS */
338 
339  // Explicit Heun-Euler.
340  case Explicit_HEUN_EULER_2_12_embedded:
341  this->alloc(2);
342  this->set_A(1, 0, 1.0);
343  this->set_B(0, 0.5);
344  this->set_B(1, 0.5);
345  this->set_B2(0, 1.0);
346  this->set_B2(1, 0.0);
347  this->set_C(0, 0.0);
348  this->set_C(1, 1.0);
349  break;
350 
351  // Explicit Bogacki-Shampine.
352  case Explicit_BOGACKI_SHAMPINE_4_23_embedded:
353  this->alloc(4);
354  this->set_A(1, 0, 1. / 2.);
355  this->set_A(3, 0, 2. / 9.);
356  this->set_A(2, 1, 3. / 4.);
357  this->set_A(3, 1, 1. / 3.);
358  this->set_A(3, 2, 4. / 9.);
359  this->set_B(0, 2. / 9.);
360  this->set_B(1, 1. / 3.);
361  this->set_B(2, 4. / 9.);
362  this->set_B(3, 0.0);
363  this->set_B2(0, 7. / 24.);
364  this->set_B2(1, 1. / 4.);
365  this->set_B2(2, 1. / 3.);
366  this->set_B2(3, 1. / 8.);
367  this->set_C(1, 1. / 2.);
368  this->set_C(2, 3. / 4.);
369  this->set_C(3, 1.0);
370  break;
371 
372  // Explicit Fehlberg.
373  case Explicit_FEHLBERG_6_45_embedded:
374  this->alloc(6);
375  this->set_A(1, 0, 1. / 4.);
376  this->set_A(2, 0, 3. / 32.);
377  this->set_A(3, 0, 1932. / 2197.);
378  this->set_A(4, 0, 439. / 216.);
379  this->set_A(5, 0, -8. / 27.);
380  this->set_A(2, 1, 9. / 32.);
381  this->set_A(3, 1, -7200. / 2197.);
382  this->set_A(4, 1, -8.);
383  this->set_A(5, 1, 2.);
384  this->set_A(3, 2, 7296. / 2197.);
385  this->set_A(4, 2, 3680. / 513.);
386  this->set_A(5, 2, -3544. / 2565.);
387  this->set_A(4, 3, -845. / 4104.);
388  this->set_A(5, 3, 1859. / 4104.);
389  this->set_A(5, 4, -11. / 40.);
390  this->set_B(0, 16. / 135.);
391  this->set_B(1, 0);
392  this->set_B(2, 6656. / 12825.);
393  this->set_B(3, 28561. / 56430.);
394  this->set_B(4, -9. / 50.);
395  this->set_B(5, 2. / 55.);
396  this->set_B2(0, 25. / 216.);
397  this->set_B2(1, 0);
398  this->set_B2(2, 1408. / 2565.);
399  this->set_B2(3, 2197. / 4104.);
400  this->set_B2(4, -1. / 5.);
401  this->set_C(1, 1. / 4.);
402  this->set_C(2, 3. / 8.);
403  this->set_C(3, 12. / 13.);
404  this->set_C(4, 1.);
405  this->set_C(5, 1. / 2.);
406  break;
407 
408  // Explicit Cash-Karp.
409  case Explicit_CASH_KARP_6_45_embedded:
410  this->alloc(6);
411  this->set_A(1, 0, 1. / 5.);
412  this->set_A(2, 0, 3. / 40.);
413  this->set_A(3, 0, 3. / 10.);
414  this->set_A(4, 0, -11. / 54.);
415  this->set_A(5, 0, 1631. / 55296.);
416  this->set_A(2, 1, 9. / 40.);
417  this->set_A(3, 1, -9. / 10.);
418  this->set_A(4, 1, 5. / 2.);
419  this->set_A(5, 1, 175. / 512.);
420  this->set_A(3, 2, 6. / 5.);
421  this->set_A(4, 2, -70. / 27.);
422  this->set_A(5, 2, 575. / 13824.);
423  this->set_A(4, 3, 35. / 27.);
424  this->set_A(5, 3, 44275. / 110592.);
425  this->set_A(5, 4, 253. / 4096.);
426  this->set_B(0, 37. / 378.);
427  this->set_B(1, 0);
428  this->set_B(2, 250. / 621.);
429  this->set_B(3, 125. / 594.);
430  this->set_B(4, 0);
431  this->set_B(5, 512. / 1771.);
432  this->set_B2(0, 2825. / 27648.);
433  this->set_B2(1, 0);
434  this->set_B2(2, 18575. / 48384.);
435  this->set_B2(3, 13525. / 55296.);
436  this->set_B2(4, 277. / 14336.);
437  this->set_B2(5, 1. / 4.);
438  this->set_C(1, 1. / 5.);
439  this->set_C(2, 3. / 10.);
440  this->set_C(3, 3. / 5.);
441  this->set_C(4, 1.);
442  this->set_C(5, 7. / 8.);
443  break;
444 
445  // Explicit Dormand-Prince.
446  case Explicit_DORMAND_PRINCE_7_45_embedded:
447  this->alloc(7);
448  this->set_A(1, 0, 1. / 5.);
449  this->set_A(2, 0, 3. / 40.);
450  this->set_A(3, 0, 44. / 45.);
451  this->set_A(4, 0, 19372. / 6561.);
452  this->set_A(5, 0, 9017. / 3168.);
453  this->set_A(6, 0, 35. / 384.);
454  this->set_A(2, 1, 9. / 40.);
455  this->set_A(3, 1, -56. / 15.);
456  this->set_A(4, 1, -25360. / 2187.);
457  this->set_A(5, 1, -355. / 33.);
458  this->set_A(6, 1, 0.);
459  this->set_A(3, 2, 32. / 9.);
460  this->set_A(4, 2, 64448. / 6561.);
461  this->set_A(5, 2, 46732. / 5247.);
462  this->set_A(6, 2, 500. / 1113.);
463  this->set_A(4, 3, -212. / 729.);
464  this->set_A(5, 3, 49. / 176.);
465  this->set_A(6, 3, 125. / 192.);
466  this->set_A(5, 4, -5103. / 18656.);
467  this->set_A(6, 4, -2187. / 6784.);
468  this->set_A(6, 5, 11. / 84.);
469  this->set_B(0, 35. / 384.);
470  this->set_B(1, 0.);
471  this->set_B(2, 500. / 1113.);
472  this->set_B(3, 125. / 192.);
473  this->set_B(4, -2187. / 6784.);
474  this->set_B(5, 11. / 84.);
475  this->set_B(6, 0.);
476  this->set_B2(0, 5179. / 57600.);
477  this->set_B2(1, 0.);
478  this->set_B2(2, 7571. / 16695.);
479  this->set_B2(3, 393. / 640.);
480  this->set_B2(4, -92097. / 339200.);
481  this->set_B2(5, 187. / 2100.);
482  this->set_B2(6, 1. / 40.);
483  this->set_C(1, 1. / 5.);
484  this->set_C(2, 3. / 10.);
485  this->set_C(3, 4. / 5.);
486  this->set_C(4, 8. / 9.);
487  this->set_C(5, 1.);
488  this->set_C(6, 1.);
489  break;
490 
491  /* EMBEDDED IMPLICIT METHODS */
492 
493  case Implicit_ESDIRK_TRBDF2_3_23_embedded:
494  this->alloc(3);
495  this->set_A(1, 0, (2 - sqrt((double)2)) / 2.);
496  this->set_A(2, 0, sqrt((double)2) / 4.);
497  this->set_A(1, 1, (2 - sqrt((double)2)) / 2.);
498  this->set_A(2, 1, sqrt((double)2) / 4.);
499  this->set_A(2, 2, (2 - sqrt((double)2)) / 2.);
500  this->set_B(0, sqrt((double)2) / 4.);
501  this->set_B(1, sqrt((double)2) / 4.);
502  this->set_B(2, (2 - sqrt((double)2)) / 2.);
503  this->set_B2(0, (1 - sqrt((double)2) / 4.) / 3.);
504  this->set_B2(1, (3 * sqrt((double)2) / 4. + 1.) / 3.);
505  this->set_B2(2, (2 - sqrt((double)2)) / 6.);
506  this->set_C(1, 2 - sqrt((double)2));
507  this->set_C(2, 1.0);
508  break;
509 
510  case Implicit_ESDIRK_TRX2_3_23_embedded:
511  this->alloc(3);
512  this->set_A(1, 0, 0.25);
513  this->set_A(2, 0, 0.25);
514  this->set_A(1, 1, 0.25);
515  this->set_A(2, 1, 0.5);
516  this->set_A(2, 2, 0.25);
517  this->set_B(0, 0.25);
518  this->set_B(1, 0.5);
519  this->set_B(2, 0.25);
520  this->set_B2(0, 1. / 6.);
521  this->set_B2(1, 2. / 3.);
522  this->set_B2(2, 1. / 6.);
523  this->set_C(1, 0.5);
524  this->set_C(2, 1.0);
525  break;
526 
527  case Implicit_SDIRK_CASH_3_23_embedded:
528  this->alloc(3);
529  this->set_A(0, 0, 0.435866521508);
530  this->set_A(1, 0, 0.2820667320);
531  this->set_A(2, 0, 1.208496649);
532  this->set_A(1, 1, 0.435866521508);
533  this->set_A(2, 1, -0.6443632015);
534  this->set_A(2, 2, 0.435866521508);
535  this->set_B(0, 1.208496649);
536  this->set_B(1, -0.6443632015);
537  this->set_B(2, 0.435866521508);
538  this->set_B2(0, 0.77263013745746);
539  this->set_B2(1, 0.22736986254254);
540  this->set_C(0, 0.435866521508);
541  this->set_C(1, 0.717933260755);
542  this->set_C(2, 1.0);
543  break;
544 
545  case Implicit_SDIRK_BILLINGTON_3_23_embedded:
546  this->alloc(3);
547  this->set_A(0, 0, 0.292893218813);
548  this->set_A(1, 0, 0.798989873223);
549  this->set_A(2, 0, 0.740789228841);
550  this->set_A(1, 1, 0.292893218813);
551  this->set_A(2, 1, 0.259210771159);
552  this->set_A(2, 2, 0.292893218813);
553  this->set_B(0, 0.691665115992);
554  this->set_B(1, 0.503597029883);
555  this->set_B(2, -0.195262145876);
556  this->set_B2(0, 0.740789228840);
557  this->set_B2(1, 0.259210771159);
558  this->set_C(0, 0.292893218813);
559  this->set_C(1, 1.091883092037);
560  this->set_C(2, 1.292893218813);
561  break;
562 
563  case Implicit_SDIRK_CASH_5_24_embedded:
564  this->alloc(5);
565  this->set_A(0, 0, 0.435866521508);
566  this->set_A(1, 0, -1.13586652150);
567  this->set_A(2, 0, 1.08543330679);
568  this->set_A(3, 0, 0.416349501547);
569  this->set_A(4, 0, 0.896869652944);
570  this->set_A(1, 1, 0.435866521508);
571  this->set_A(2, 1, -0.721299828287);
572  this->set_A(3, 1, 0.190984004184);
573  this->set_A(4, 1, 0.0182725272734);
574  this->set_A(2, 2, 0.435866521508);
575  this->set_A(3, 2, -0.118643265417);
576  this->set_A(4, 2, -0.0845900310706);
577  this->set_A(3, 3, 0.435866521508);
578  this->set_A(4, 3, -0.266418670647);
579  this->set_A(4, 4, 0.435866521508);
580  this->set_B(0, 0.896869652944);
581  this->set_B(1, 0.0182725272734);
582  this->set_B(2, -0.0845900310706);
583  this->set_B(3, -0.266418670647);
584  this->set_B(4, 0.435866521508);
585  this->set_B2(0, (-0.7 - 0.5) / (-0.7 - 0.435866521508));
586  this->set_B2(1, (0.5 - 0.435866521508) / (-0.7 - 0.435866521508));
587  this->set_B2(2, 0.0);
588  this->set_B2(3, 0.0);
589  this->set_B2(4, 0.0);
590  this->set_C(0, 0.435866521508);
591  this->set_C(1, -0.7);
592  this->set_C(2, 0.8);
593  this->set_C(3, 0.924556761814);
594  this->set_C(4, 1.0);
595  break;
596 
597  case Implicit_SDIRK_CASH_5_34_embedded:
598  this->alloc(5);
599  this->set_A(0, 0, 0.435866521508);
600  this->set_A(1, 0, -1.13586652150);
601  this->set_A(2, 0, 1.08543330679);
602  this->set_A(3, 0, 0.416349501547);
603  this->set_A(4, 0, 0.896869652944);
604  this->set_A(1, 1, 0.435866521508);
605  this->set_A(2, 1, -0.721299828287);
606  this->set_A(3, 1, 0.190984004184);
607  this->set_A(4, 1, 0.0182725272734);
608  this->set_A(2, 2, 0.435866521508);
609  this->set_A(3, 2, -0.118643265417);
610  this->set_A(4, 2, -0.0845900310706);
611  this->set_A(3, 3, 0.435866521508);
612  this->set_A(4, 3, -0.266418670647);
613  this->set_A(4, 4, 0.435866521508);
614  this->set_B(0, 0.896869652944);
615  this->set_B(1, 0.0182725272734);
616  this->set_B(2, -0.0845900310706);
617  this->set_B(3, -0.266418670647);
618  this->set_B(4, 0.435866521508);
619  this->set_B2(0, 0.776691932910);
620  this->set_B2(1, 0.0297472791484);
621  this->set_B2(2, -0.0267440239074);
622  this->set_B2(3, 0.220304811849);
623  this->set_B2(4, 0.0);
624  this->set_C(0, 0.435866521508);
625  this->set_C(1, -0.7);
626  this->set_C(2, 0.8);
627  this->set_C(3, 0.924556761814);
628  this->set_C(4, 1.0);
629  break;
630 
631  case Implicit_DIRK_ISMAIL_7_45_embedded: // Implicit embedded DIRK method with orders 4 and 5.
632  this->alloc(7);
633  this->set_A(0, 0, 0);
634  this->set_A(1, 0, 0.28589);
635  this->set_A(2, 0, 0.142945);
636  this->set_A(3, 0, 0.16803599);
637  this->set_A(4, 0, 0.182315);
638  this->set_A(5, 0, 0.24756392);
639  this->set_A(6, 0, 0.13001804);
640  this->set_A(1, 1, 0.28589);
641  this->set_A(2, 1, 0.924011005);
642  this->set_A(3, 1, -0.049416510);
643  this->set_A(4, 1, -0.112951603);
644  this->set_A(5, 1, -0.425378071);
645  this->set_A(6, 1, 0);
646  this->set_A(2, 2, 0.28589);
647  this->set_A(3, 2, -0.004509476);
648  this->set_A(4, 2, -0.027793233);
649  this->set_A(5, 2, -0.107036282);
650  this->set_A(6, 2, -0.019290177);
651  this->set_A(3, 3, 0.28589);
652  this->set_A(4, 3, 0.422539833);
653  this->set_A(5, 3, 0.395700134);
654  this->set_A(6, 3, 0.535386266);
655  this->set_A(4, 4, 0.28589);
656  this->set_A(5, 4, 0.503260302);
657  this->set_A(6, 4, 0.234313169);
658  this->set_A(5, 5, 0.28589);
659  this->set_A(6, 5, -0.166317293);
660  this->set_A(6, 6, 0.28589);
661  this->set_B(0, 0.13001804);
662  this->set_B(1, 0);
663  this->set_B(2, -0.019290177);
664  this->set_B(3, 0.535386266);
665  this->set_B(4, 0.234313169);
666  this->set_B(5, -0.166317293);
667  this->set_B(6, 0.28589);
668  this->set_B2(0, 0.094388663);
669  this->set_B2(1, 0);
670  this->set_B2(2, -0.039782614);
671  this->set_B2(3, 0.745608552);
672  this->set_B2(4, -0.505129807);
673  this->set_B2(5, 0.704915206);
674  this->set_B2(6, 0);
675  this->set_C(0, 0);
676  this->set_C(1, 0.57178);
677  this->set_C(2, 1.352846);
678  this->set_C(3, 0.4);
679  this->set_C(4, 0.75);
680  this->set_C(5, 0.9);
681  this->set_C(6, 1.0);
682  break;
683 
684  default: throw Hermes::Exceptions::Exception("Unknown Butcher's table.");
685  }
686  }
687 
688  void ButcherTable::alloc(unsigned int size)
689  {
690  // Size.
691  this->size = size;
692  // A array.
693  this->A = new_matrix<double>(size, size);
694  for (unsigned int i = 0; i < size; i++)
695  {
696  for (unsigned int j = 0; j < size; j++) this->A[i][j] = 0;
697  }
698  // B array.
699  this->B = malloc_with_check<ButcherTable, double>(size, this);
700  for (unsigned int j = 0; j < size; j++) this->B[j] = 0;
701  // B2 array.
702  this->B2 = malloc_with_check<ButcherTable, double>(size, this);
703  for (unsigned int j = 0; j < size; j++) this->B2[j] = 0;
704  // C array.
705  this->C = malloc_with_check<ButcherTable, double>(size, this);
706  for (unsigned int j = 0; j < size; j++) this->C[j] = 0;
707  }
708 
709  ButcherTable::~ButcherTable()
710  {
711  this->free();
712  }
713 
714  void ButcherTable::free()
715  {
716  free_with_check(this->B);
717  free_with_check(this->B2);
718  free_with_check(this->C);
719  }
720 
721  double ButcherTable::get_B(unsigned int i)
722  {
723  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
724  return this->B[i];
725  }
726 
727  double ButcherTable::get_B2(unsigned int i)
728  {
729  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
730  return this->B2[i];
731  }
732 
733  double ButcherTable::get_C(unsigned int i)
734  {
735  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
736  return this->C[i];
737  }
738 
739  void ButcherTable::set_B(unsigned int i, double val)
740  {
741  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
742  this->B[i] = val;
743  }
744 
745  void ButcherTable::set_B2(unsigned int i, double val)
746  {
747  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
748  this->B2[i] = val;
749  }
750 
751  void ButcherTable::set_C(unsigned int i, double val)
752  {
753  if (i > size) throw Hermes::Exceptions::Exception("Invalid access to a Butcher's table.");
754  this->C[i] = val;
755  }
756 
757  bool ButcherTable::is_explicit()
758  {
759  bool result = true;
760  for (unsigned int i = 0; i < size; i++)
761  {
762  for (unsigned int j = 0; j < size; j++)
763  {
764  double val_ij = get_A(i, j);
765  if (j >= i && fabs(val_ij) > Hermes::HermesSqrtEpsilon) result = false;
766  }
767  }
768 
769  return result;
770  }
771 
772  bool ButcherTable::is_diagonally_implicit()
773  {
774  bool result = true;
775  for (unsigned int i = 0; i < size; i++)
776  {
777  for (unsigned int j = 0; j < size; j++)
778  {
779  double val_ij = get_A(i, j);
780  if (j > i && fabs(val_ij) > Hermes::HermesSqrtEpsilon) result = false;
781  }
782  }
783 
784  return result;
785  }
786 
787  bool ButcherTable::is_fully_implicit()
788  {
789  bool result = false;
790  for (unsigned int i = 0; i < size; i++)
791  {
792  for (unsigned int j = 0; j < size; j++)
793  {
794  double val_ij = get_A(i, j);
795  if (j > i && fabs(val_ij) > Hermes::HermesSqrtEpsilon) result = true;
796  }
797  }
798 
799  return result;
800  }
801 
802  bool ButcherTable::is_embedded()
803  {
804  // Test whether B2 row is not zero.
805  double sum = 0;
806  for (unsigned int i = 0; i < size; i++) sum += fabs(B2[i]);
807  if (sum < Hermes::HermesSqrtEpsilon) return false;
808  else return true;
809  }
810 
811  void ButcherTable::switch_B_rows()
812  {
813  // Test whether nonzero B2 row exists.
814  if (this->is_embedded() == false)
815  throw Hermes::Exceptions::Exception("ButcherTable::switch_B_rows(): Zero B2 row detected.");
816 
817  // Switch B rows.
818  for (unsigned int i = 0; i < size; i++)
819  {
820  double tmp = B[i];
821  B[i] = B2[i];
822  B2[i] = tmp;
823  }
824  }
825 }
Hermes
General namespace for the Hermes library.
Definition: algebra_mixins.h:30
Hermes::Implicit_Lobatto_IIIC_2_2
Implicit Lobatto IIIB-2 method.
Definition: table.h:50
Hermes::Exceptions::Exception
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
Hermes::Implicit_RK_1
Implicit Runge-Kutta RK-1 (implicit Euler).
Definition: table.h:43
Hermes::Implicit_SIRK_2_2
Implicit SIRK-2-2 method.
Definition: table.h:45
Hermes::ButcherTable::switch_B_rows
void switch_B_rows()
Definition: table.cpp:811
Hermes::Implicit_Lobatto_IIIA_3_4
Implicit Lobatto IIIA-4 method.
Definition: table.h:51
Hermes::Explicit_HEUN_EULER_2_12_embedded
Heun-Euler (orders 1 and 2),.
Definition: table.h:59
Hermes::Explicit_FEHLBERG_6_45_embedded
Fehlberg (orders 4 and 5),.
Definition: table.h:61
Hermes::Explicit_BOGACKI_SHAMPINE_4_23_embedded
Bogacki-Shampine (orders 2 and 3),.
Definition: table.h:60
memory_handling.h
File containing common definitions, and basic global enums etc. for HermesCommon. ...
Hermes::Explicit_RK_1
Explicit Runge-Kutta RK-1 (explicit Euler).
Definition: table.h:36
Hermes::Implicit_ESDIRK_TRBDF2_3_23_embedded
of TR-BDF2 by Hosea and Shampine.
Definition: table.h:67
Hermes::Implicit_SDIRK_5_4
Implicit SDIRK-5-4 method.
Definition: table.h:55
Hermes::Explicit_CASH_KARP_6_45_embedded
Cash-Karp (orders 4 and 5),.
Definition: table.h:62
Hermes::Implicit_Lobatto_IIIA_2_2
Implicit Lobatto IIIA-2 method.
Definition: table.h:48
Hermes::Explicit_RK_3
Explicit Runge-Kutta RK-3 method.
Definition: table.h:38
Hermes::Implicit_Radau_IIA_3_5
Implicit Radau IIA-5 method.
Definition: table.h:54
Hermes::Algebra::DenseMatrixOperations
Contains operation on dense matrices.
Definition: dense_matrix_operations.h:36
Hermes::Implicit_SDIRK_CASH_3_23_embedded
with Error Estimates by J.R. Cash
Definition: table.h:71
Hermes::Explicit_DORMAND_PRINCE_7_45_embedded
Dormand-Prince (orders 4 and 5),.
Definition: table.h:63
Hermes::Explicit_RK_4
Explicit Runge-Kutta RK-4 method.
Definition: table.h:39
Hermes::Implicit_ESDIRK_TRX2_3_23_embedded
of TR-BDF2 by Hosea and Shampine.
Definition: table.h:69
Hermes::ButcherTable::C
double * C
on embedded R-K methods.
Definition: table.h:130
Hermes::ButcherTable::B2
double * B2
This is the second B-row for adaptivity based.
Definition: table.h:128
matrix.h
Basic matrix classes and operations.
Hermes::Explicit_RK_2
Explicit Runge-Kutta RK-2 method.
Definition: table.h:37
Hermes::ButcherTableType
ButcherTableType
Butcher's tables type. The last number in the name always means order, the one before last (if provid...
Definition: table.h:32
Hermes::Implicit_SDIRK_2_2
Implicit SDIRK-2-2 method.
Definition: table.h:47
Hermes::Implicit_SDIRK_CASH_5_34_embedded
with Error Estimates by J.R. Cash
Definition: table.h:77
Hermes::Implicit_Lobatto_IIIB_2_2
Implicit Lobatto IIIB-2 method.
Definition: table.h:49
Hermes::Implicit_Lobatto_IIIC_3_4
Implicit Lobatto IIIB-4 method.
Definition: table.h:53
dense_matrix_operations.h
Dense (small) simply stored matrix operations.
Hermes::Implicit_SDIRK_CASH_5_24_embedded
with Error Estimates by J.R. Cash
Definition: table.h:75
Hermes::Implicit_SDIRK_BILLINGTON_3_23_embedded
Codes for the Solution of Stiff and Non-Stiff systems of ODE's.
Definition: table.h:73
Hermes::Implicit_ESIRK_2_2
Implicit ESIRK-2-2 method.
Definition: table.h:46
Hermes::Implicit_Crank_Nicolson_2_2
Implicit Crank_Nicolson method.
Definition: table.h:44
Hermes::Implicit_Lobatto_IIIB_3_4
Implicit Lobatto IIIB-4 method.
Definition: table.h:52