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