Rythmos - Transient Integration for Differential Equations  Version of the Day
Rythmos_RKButcherTableau.hpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Rythmos Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25 //
26 // ***********************************************************************
27 //@HEADER
28 
29 
30 #ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HPP
31 #define RYTHMOS_RK_BUTCHER_TABLEAU_HPP
32 
33 // disable clang warnings
34 #ifdef __clang__
35 #pragma clang system_header
36 #endif
37 
38 #include "Rythmos_Types.hpp"
39 #include "Rythmos_RKButcherTableauBase.hpp"
40 
41 #include "Teuchos_Assert.hpp"
42 #include "Teuchos_as.hpp"
43 #include "Teuchos_StandardParameterEntryValidators.hpp"
44 #include "Teuchos_Describable.hpp"
45 #include "Teuchos_VerboseObject.hpp"
46 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
47 #include "Teuchos_ParameterListAcceptor.hpp"
48 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
49 
50 #include "Thyra_ProductVectorBase.hpp"
51 
52 namespace Rythmos {
53 
54  using Teuchos::as;
55 
56  inline const std::string RKBT_ForwardEuler_name() { return "Forward Euler"; } // done
57  inline const std::string RKBT_BackwardEuler_name() { return "Backward Euler"; } // done
58  inline const std::string Explicit4Stage_name() { return "Explicit 4 Stage"; } // done
59  inline const std::string Explicit3_8Rule_name() { return "Explicit 3/8 Rule"; } // done
60 
61  inline const std::string ExplicitTrapezoidal_name() { return "Explicit Trapezoidal"; } // done
62  inline const std::string Explicit2Stage2ndOrderRunge_name() { return "Explicit 2 Stage 2nd order by Runge"; } // done
63  inline const std::string Explicit3Stage3rdOrderHeun_name() { return "Explicit 3 Stage 3rd order by Heun"; } // done
64  inline const std::string Explicit3Stage3rdOrder_name() { return "Explicit 3 Stage 3rd order"; } // done
65  inline const std::string Explicit3Stage3rdOrderTVD_name() { return "Explicit 3 Stage 3rd order TVD"; } // done
66  inline const std::string Explicit4Stage3rdOrderRunge_name() { return "Explicit 4 Stage 3rd order by Runge"; } // done
67 
68  inline const std::string IRK1StageTheta_name() { return "IRK 1 Stage Theta Method"; } // done
69  inline const std::string IRK2StageTheta_name() { return "IRK 2 Stage Theta Method"; } // done
70  inline const std::string Implicit1Stage2ndOrderGauss_name() { return "Implicit 1 Stage 2nd order Gauss"; } // done
71  inline const std::string Implicit2Stage4thOrderGauss_name() { return "Implicit 2 Stage 4th order Gauss"; } // done
72  inline const std::string Implicit3Stage6thOrderGauss_name() { return "Implicit 3 Stage 6th order Gauss"; } // done
73 
74  inline const std::string Implicit1Stage1stOrderRadauA_name() { return "Implicit 1 Stage 1st order Radau left"; } // done
75  inline const std::string Implicit2Stage3rdOrderRadauA_name() { return "Implicit 2 Stage 3rd order Radau left"; } // done
76  inline const std::string Implicit3Stage5thOrderRadauA_name() { return "Implicit 3 Stage 5th order Radau left"; } // done
77 
78  inline const std::string Implicit1Stage1stOrderRadauB_name() { return "Implicit 1 Stage 1st order Radau right"; } // done
79  inline const std::string Implicit2Stage3rdOrderRadauB_name() { return "Implicit 2 Stage 3rd order Radau right"; } // done
80  inline const std::string Implicit3Stage5thOrderRadauB_name() { return "Implicit 3 Stage 5th order Radau right"; } // done
81 
82  inline const std::string Implicit2Stage2ndOrderLobattoA_name() { return "Implicit 2 Stage 2nd order Lobatto A"; } // done
83  inline const std::string Implicit3Stage4thOrderLobattoA_name() { return "Implicit 3 Stage 4th order Lobatto A"; } // done
84  inline const std::string Implicit4Stage6thOrderLobattoA_name() { return "Implicit 4 Stage 6th order Lobatto A"; } // done
85 
86  inline const std::string Implicit2Stage2ndOrderLobattoB_name() { return "Implicit 2 Stage 2nd order Lobatto B"; } // done
87  inline const std::string Implicit3Stage4thOrderLobattoB_name() { return "Implicit 3 Stage 4th order Lobatto B"; } // done
88  inline const std::string Implicit4Stage6thOrderLobattoB_name() { return "Implicit 4 Stage 6th order Lobatto B"; } // done
89 
90  inline const std::string Implicit2Stage2ndOrderLobattoC_name() { return "Implicit 2 Stage 2nd order Lobatto C"; } // done
91  inline const std::string Implicit3Stage4thOrderLobattoC_name() { return "Implicit 3 Stage 4th order Lobatto C"; } // done
92  inline const std::string Implicit4Stage6thOrderLobattoC_name() { return "Implicit 4 Stage 6th order Lobatto C"; } // done
93 
94  inline const std::string Implicit2Stage4thOrderHammerHollingsworth_name() { return "Implicit 2 Stage 4th Order Hammer & Hollingsworth"; } // done
95  inline const std::string Implicit3Stage6thOrderKuntzmannButcher_name() { return "Implicit 3 Stage 6th Order Kuntzmann & Butcher"; } // done
96  inline const std::string Implicit4Stage8thOrderKuntzmannButcher_name() { return "Implicit 4 Stage 8th Order Kuntzmann & Butcher"; } // done
97 
98  inline const std::string DIRK2Stage3rdOrder_name() { return "Diagonal IRK 2 Stage 3rd order"; } // done
99 
100  inline const std::string SDIRK2Stage2ndOrder_name() { return "Singly Diagonal IRK 2 Stage 2nd order"; } // done
101  inline const std::string SDIRK2Stage3rdOrder_name() { return "Singly Diagonal IRK 2 Stage 3rd order"; } // done
102  inline const std::string SDIRK5Stage5thOrder_name() { return "Singly Diagonal IRK 5 Stage 5th order"; } // done
103  inline const std::string SDIRK5Stage4thOrder_name() { return "Singly Diagonal IRK 5 Stage 4th order"; } // done
104  inline const std::string SDIRK3Stage4thOrder_name() { return "Singly Diagonal IRK 3 Stage 4th order"; } // done
105 
106 template<class Scalar>
107 class RKButcherTableauDefaultBase :
108  virtual public RKButcherTableauBase<Scalar>,
109  virtual public Teuchos::ParameterListAcceptorDefaultBase
110 {
111  public:
113  virtual int numStages() const { return A_.numRows(); }
115  virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const { return A_; }
117  virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const { return b_; }
119  virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const { return c_; }
121  virtual int order() const { return order_; }
123  virtual void setDescription(std::string longDescription) { longDescription_ = longDescription; }
124 
126  virtual void initialize(
127  const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
128  const Teuchos::SerialDenseVector<int,Scalar>& b_in,
129  const Teuchos::SerialDenseVector<int,Scalar>& c_in,
130  const int order_in,
131  const std::string& longDescription_in
132  )
133  {
134  const int numStages_in = A_in.numRows();
135  TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
136  TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
137  TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
138  TEUCHOS_ASSERT_EQUALITY( c_in.length(), numStages_in );
139  TEUCHOS_ASSERT( order_in > 0 );
140  A_ = A_in;
141  b_ = b_in;
142  c_ = c_in;
143  order_ = order_in;
144  longDescription_ = longDescription_in;
145  }
146 
147  /* \brief Redefined from Teuchos::ParameterListAcceptorDefaultBase */
149 
151  virtual void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
152  {
153  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
154  paramList->validateParameters(*this->getValidParameters());
155  Teuchos::readVerboseObjectSublist(&*paramList,this);
156  setMyParamList(paramList);
157  }
158 
160  virtual RCP<const Teuchos::ParameterList> getValidParameters() const
161  {
162  if (is_null(validPL_)) {
163  validPL_ = Teuchos::parameterList();
164  validPL_->set("Description","",this->getMyDescription());
165  Teuchos::setupVerboseObjectSublist(&*validPL_);
166  }
167  return validPL_;
168  }
169 
171 
172  /* \brief Redefined from Teuchos::Describable */
174 
176  virtual std::string description() const { return "Rythmos::RKButcherTableauDefaultBase"; }
177 
179  virtual void describe(
180  Teuchos::FancyOStream &out,
181  const Teuchos::EVerbosityLevel verbLevel
182  ) const
183  {
184  if (verbLevel != Teuchos::VERB_NONE) {
185  out << this->description() << std::endl;
186  out << this->getMyDescription() << std::endl;
187  out << "number of Stages = " << this->numStages() << std::endl;
188  out << "A = " << this->A() << std::endl;
189  out << "b = " << this->b() << std::endl;
190  out << "c = " << this->c() << std::endl;
191  out << "order = " << this->order() << std::endl;
192  }
193  }
194 
196 
197  protected:
198  void setMyDescription(std::string longDescription) { longDescription_ = longDescription; }
199  const std::string& getMyDescription() const { return longDescription_; }
200 
201  void setMy_A(const Teuchos::SerialDenseMatrix<int,Scalar>& new_A) { A_ = new_A; }
202  void setMy_b(const Teuchos::SerialDenseVector<int,Scalar>& new_b) { b_ = new_b; }
203  void setMy_c(const Teuchos::SerialDenseVector<int,Scalar>& new_c) { c_ = new_c; }
204  void setMy_order(const int& new_order) { order_ = new_order; }
205 
206  void setMyValidParameterList( const RCP<ParameterList> validPL ) { validPL_ = validPL; }
207  RCP<ParameterList> getMyNonconstValidParameterList() { return validPL_; }
208 
209  private:
210  Teuchos::SerialDenseMatrix<int,Scalar> A_;
211  Teuchos::SerialDenseVector<int,Scalar> b_;
212  Teuchos::SerialDenseVector<int,Scalar> c_;
213  int order_;
214  std::string longDescription_;
215  mutable RCP<ParameterList> validPL_;
216 };
217 
218 
219 // Nonmember constructor
220 template<class Scalar>
221 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
222 {
223  return(rcp(new RKButcherTableauDefaultBase<Scalar>()));
224 }
225 
226 // Nonmember constructor
227 template<class Scalar>
228 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau(
229  const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
230  const Teuchos::SerialDenseVector<int,Scalar>& b_in,
231  const Teuchos::SerialDenseVector<int,Scalar>& c_in,
232  int order_in,
233  const std::string& description_in = ""
234  )
235 {
236  RCP<RKButcherTableauDefaultBase<Scalar> > rkbt = rcp(new RKButcherTableauDefaultBase<Scalar>());
237  rkbt->initialize(A_in,b_in,c_in,order_in,description_in);
238  return(rkbt);
239 }
240 
241 
242 template<class Scalar>
243 class BackwardEuler_RKBT :
244  virtual public RKButcherTableauDefaultBase<Scalar>
245 {
246  public:
247  BackwardEuler_RKBT()
248  {
249  std::ostringstream myDescription;
250  myDescription << RKBT_BackwardEuler_name() << "\n"
251  << "c = [ 1 ]'\n"
252  << "A = [ 1 ]\n"
253  << "b = [ 1 ]'" << std::endl;
254  typedef ScalarTraits<Scalar> ST;
255  Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
256  myA(0,0) = ST::one();
257  Teuchos::SerialDenseVector<int,Scalar> myb(1);
258  myb(0) = ST::one();
259  Teuchos::SerialDenseVector<int,Scalar> myc(1);
260  myc(0) = ST::one();
261 
262  this->setMyDescription(myDescription.str());
263  this->setMy_A(myA);
264  this->setMy_b(myb);
265  this->setMy_c(myc);
266  this->setMy_order(1);
267  }
268 };
269 
270 
271 
272 template<class Scalar>
273 class ForwardEuler_RKBT :
274  virtual public RKButcherTableauDefaultBase<Scalar>
275 {
276  public:
277 
278  ForwardEuler_RKBT()
279  {
280  std::ostringstream myDescription;
281  myDescription << RKBT_ForwardEuler_name() << "\n"
282  << "c = [ 0 ]'\n"
283  << "A = [ 0 ]\n"
284  << "b = [ 1 ]'" << std::endl;
285  typedef ScalarTraits<Scalar> ST;
286  Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
287  Teuchos::SerialDenseVector<int,Scalar> myb(1);
288  myb(0) = ST::one();
289  Teuchos::SerialDenseVector<int,Scalar> myc(1);
290 
291  this->setMyDescription(myDescription.str());
292  this->setMy_A(myA);
293  this->setMy_b(myb);
294  this->setMy_c(myc);
295  this->setMy_order(1);
296  }
297 };
298 
299 
300 template<class Scalar>
301 class Explicit4Stage4thOrder_RKBT :
302  virtual public RKButcherTableauDefaultBase<Scalar>
303 {
304  public:
305  Explicit4Stage4thOrder_RKBT()
306  {
307  std::ostringstream myDescription;
308  myDescription << Explicit4Stage_name() << "\n"
309  << "\"The\" Runge-Kutta Method (explicit):\n"
310  << "Solving Ordinary Differential Equations I:\n"
311  << "Nonstiff Problems, 2nd Revised Edition\n"
312  << "E. Hairer, S.P. Norsett, G. Wanner\n"
313  << "Table 1.2, pg 138\n"
314  << "c = [ 0 1/2 1/2 1 ]'\n"
315  << "A = [ 0 ] \n"
316  << " [ 1/2 0 ]\n"
317  << " [ 0 1/2 0 ]\n"
318  << " [ 0 0 1 0 ]\n"
319  << "b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
320  typedef ScalarTraits<Scalar> ST;
321  Scalar one = ST::one();
322  Scalar zero = ST::zero();
323  Scalar onehalf = ST::one()/(2*ST::one());
324  Scalar onesixth = ST::one()/(6*ST::one());
325  Scalar onethird = ST::one()/(3*ST::one());
326 
327  int myNumStages = 4;
328  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
329  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
330  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
331 
332  // Fill A:
333  myA(0,0) = zero;
334  myA(0,1) = zero;
335  myA(0,2) = zero;
336  myA(0,3) = zero;
337 
338  myA(1,0) = onehalf;
339  myA(1,1) = zero;
340  myA(1,2) = zero;
341  myA(1,3) = zero;
342 
343  myA(2,0) = zero;
344  myA(2,1) = onehalf;
345  myA(2,2) = zero;
346  myA(2,3) = zero;
347 
348  myA(3,0) = zero;
349  myA(3,1) = zero;
350  myA(3,2) = one;
351  myA(3,3) = zero;
352 
353  // Fill myb:
354  myb(0) = onesixth;
355  myb(1) = onethird;
356  myb(2) = onethird;
357  myb(3) = onesixth;
358 
359  // fill b_c_
360  myc(0) = zero;
361  myc(1) = onehalf;
362  myc(2) = onehalf;
363  myc(3) = one;
364 
365  this->setMyDescription(myDescription.str());
366  this->setMy_A(myA);
367  this->setMy_b(myb);
368  this->setMy_c(myc);
369  this->setMy_order(4);
370  }
371 };
372 
373 
374 template<class Scalar>
375 class Explicit3_8Rule_RKBT :
376  virtual public RKButcherTableauDefaultBase<Scalar>
377 {
378  public:
379  Explicit3_8Rule_RKBT()
380  {
381 
382  std::ostringstream myDescription;
383  myDescription << Explicit3_8Rule_name() << "\n"
384  << "Solving Ordinary Differential Equations I:\n"
385  << "Nonstiff Problems, 2nd Revised Edition\n"
386  << "E. Hairer, S.P. Norsett, G. Wanner\n"
387  << "Table 1.2, pg 138\n"
388  << "c = [ 0 1/3 2/3 1 ]'\n"
389  << "A = [ 0 ]\n"
390  << " [ 1/3 0 ]\n"
391  << " [-1/3 1 0 ]\n"
392  << " [ 1 -1 1 0 ]\n"
393  << "b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
394  typedef ScalarTraits<Scalar> ST;
395  int myNumStages = 4;
396  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
397  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
398  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
399 
400  Scalar one = ST::one();
401  Scalar zero = ST::zero();
402  Scalar one_third = as<Scalar>(ST::one()/(3*ST::one()));
403  Scalar two_third = as<Scalar>(2*ST::one()/(3*ST::one()));
404  Scalar one_eighth = as<Scalar>(ST::one()/(8*ST::one()));
405  Scalar three_eighth = as<Scalar>(3*ST::one()/(8*ST::one()));
406 
407  // Fill myA:
408  myA(0,0) = zero;
409  myA(0,1) = zero;
410  myA(0,2) = zero;
411  myA(0,3) = zero;
412 
413  myA(1,0) = one_third;
414  myA(1,1) = zero;
415  myA(1,2) = zero;
416  myA(1,3) = zero;
417 
418  myA(2,0) = as<Scalar>(-one_third);
419  myA(2,1) = one;
420  myA(2,2) = zero;
421  myA(2,3) = zero;
422 
423  myA(3,0) = one;
424  myA(3,1) = as<Scalar>(-one);
425  myA(3,2) = one;
426  myA(3,3) = zero;
427 
428  // Fill myb:
429  myb(0) = one_eighth;
430  myb(1) = three_eighth;
431  myb(2) = three_eighth;
432  myb(3) = one_eighth;
433 
434  // Fill myc:
435  myc(0) = zero;
436  myc(1) = one_third;
437  myc(2) = two_third;
438  myc(3) = one;
439 
440  this->setMyDescription(myDescription.str());
441  this->setMy_A(myA);
442  this->setMy_b(myb);
443  this->setMy_c(myc);
444  this->setMy_order(4);
445  }
446 };
447 
448 
449 template<class Scalar>
450 class Explicit4Stage3rdOrderRunge_RKBT :
451  virtual public RKButcherTableauDefaultBase<Scalar>
452 {
453  public:
454  Explicit4Stage3rdOrderRunge_RKBT()
455  {
456 
457  std::ostringstream myDescription;
458  myDescription << Explicit4Stage3rdOrderRunge_name() << "\n"
459  << "Solving Ordinary Differential Equations I:\n"
460  << "Nonstiff Problems, 2nd Revised Edition\n"
461  << "E. Hairer, S.P. Norsett, G. Wanner\n"
462  << "Table 1.1, pg 135\n"
463  << "c = [ 0 1/2 1 1 ]'\n"
464  << "A = [ 0 ]\n"
465  << " [ 1/2 0 ]\n"
466  << " [ 0 1 0 ]\n"
467  << " [ 0 0 1 0 ]\n"
468  << "b = [ 1/6 2/3 0 1/6 ]'" << std::endl;
469  typedef ScalarTraits<Scalar> ST;
470  int myNumStages = 4;
471  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
472  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
473  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
474 
475  Scalar one = ST::one();
476  Scalar onehalf = ST::one()/(2*ST::one());
477  Scalar onesixth = ST::one()/(6*ST::one());
478  Scalar twothirds = 2*ST::one()/(3*ST::one());
479  Scalar zero = ST::zero();
480 
481  // Fill A:
482  myA(0,0) = zero;
483  myA(0,1) = zero;
484  myA(0,2) = zero;
485  myA(0,3) = zero;
486 
487  myA(1,0) = onehalf;
488  myA(1,1) = zero;
489  myA(1,2) = zero;
490  myA(1,3) = zero;
491 
492  myA(2,0) = zero;
493  myA(2,1) = one;
494  myA(2,2) = zero;
495  myA(2,3) = zero;
496 
497  myA(3,0) = zero;
498  myA(3,1) = zero;
499  myA(3,2) = one;
500  myA(3,3) = zero;
501 
502  // Fill b:
503  myb(0) = onesixth;
504  myb(1) = twothirds;
505  myb(2) = zero;
506  myb(3) = onesixth;
507 
508  // Fill myc:
509  myc(0) = zero;
510  myc(1) = onehalf;
511  myc(2) = one;
512  myc(3) = one;
513 
514  this->setMyDescription(myDescription.str());
515  this->setMy_A(myA);
516  this->setMy_b(myb);
517  this->setMy_c(myc);
518  this->setMy_order(3);
519  }
520 };
521 
522 
523 template<class Scalar>
524 class Explicit3Stage3rdOrder_RKBT :
525  virtual public RKButcherTableauDefaultBase<Scalar>
526 {
527  public:
528  Explicit3Stage3rdOrder_RKBT()
529  {
530 
531  std::ostringstream myDescription;
532  myDescription << Explicit3Stage3rdOrder_name() << "\n"
533  << "c = [ 0 1/2 1 ]'\n"
534  << "A = [ 0 ]\n"
535  << " [ 1/2 0 ]\n"
536  << " [ -1 2 0 ]\n"
537  << "b = [ 1/6 4/6 1/6 ]'" << std::endl;
538  typedef ScalarTraits<Scalar> ST;
539  Scalar one = ST::one();
540  Scalar two = as<Scalar>(2*ST::one());
541  Scalar zero = ST::zero();
542  Scalar onehalf = ST::one()/(2*ST::one());
543  Scalar onesixth = ST::one()/(6*ST::one());
544  Scalar foursixth = 4*ST::one()/(6*ST::one());
545 
546  int myNumStages = 3;
547  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
548  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
549  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
550 
551  // Fill myA:
552  myA(0,0) = zero;
553  myA(0,1) = zero;
554  myA(0,2) = zero;
555 
556  myA(1,0) = onehalf;
557  myA(1,1) = zero;
558  myA(1,2) = zero;
559 
560  myA(2,0) = -one;
561  myA(2,1) = two;
562  myA(2,2) = zero;
563 
564  // Fill myb:
565  myb(0) = onesixth;
566  myb(1) = foursixth;
567  myb(2) = onesixth;
568 
569  // fill b_c_
570  myc(0) = zero;
571  myc(1) = onehalf;
572  myc(2) = one;
573 
574  this->setMyDescription(myDescription.str());
575  this->setMy_A(myA);
576  this->setMy_b(myb);
577  this->setMy_c(myc);
578  this->setMy_order(3);
579  }
580 };
581 
582 
583 template<class Scalar>
584 class Explicit3Stage3rdOrderTVD_RKBT :
585  virtual public RKButcherTableauDefaultBase<Scalar>
586 {
587  public:
588  Explicit3Stage3rdOrderTVD_RKBT()
589  {
590 
591  std::ostringstream myDescription;
592  myDescription << Explicit3Stage3rdOrderTVD_name() << "\n"
593  << "Sigal Gottlieb and Chi-Wang Shu\n"
594  << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
595  << "Mathematics of Computation\n"
596  << "Volume 67, Number 221, January 1998, pp. 73-85\n"
597  << "c = [ 0 1 1/2 ]'\n"
598  << "A = [ 0 ]\n"
599  << " [ 1 0 ]\n"
600  << " [ 1/4 1/4 0 ]\n"
601  << "b = [ 1/6 1/6 4/6 ]'\n"
602  << "This is also written in the following set of updates.\n"
603  << "u1 = u^n + dt L(u^n)\n"
604  << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
605  << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3"
606  << std::endl;
607  typedef ScalarTraits<Scalar> ST;
608  Scalar one = ST::one();
609  Scalar zero = ST::zero();
610  Scalar onehalf = ST::one()/(2*ST::one());
611  Scalar onefourth = ST::one()/(4*ST::one());
612  Scalar onesixth = ST::one()/(6*ST::one());
613  Scalar foursixth = 4*ST::one()/(6*ST::one());
614 
615  int myNumStages = 3;
616  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
617  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
618  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
619 
620  // Fill myA:
621  myA(0,0) = zero;
622  myA(0,1) = zero;
623  myA(0,2) = zero;
624 
625  myA(1,0) = one;
626  myA(1,1) = zero;
627  myA(1,2) = zero;
628 
629  myA(2,0) = onefourth;
630  myA(2,1) = onefourth;
631  myA(2,2) = zero;
632 
633  // Fill myb:
634  myb(0) = onesixth;
635  myb(1) = onesixth;
636  myb(2) = foursixth;
637 
638  // fill b_c_
639  myc(0) = zero;
640  myc(1) = one;
641  myc(2) = onehalf;
642 
643  this->setMyDescription(myDescription.str());
644  this->setMy_A(myA);
645  this->setMy_b(myb);
646  this->setMy_c(myc);
647  this->setMy_order(3);
648  }
649 };
650 
651 
652 template<class Scalar>
653 class Explicit3Stage3rdOrderHeun_RKBT :
654  virtual public RKButcherTableauDefaultBase<Scalar>
655 {
656  public:
657  Explicit3Stage3rdOrderHeun_RKBT()
658  {
659  std::ostringstream myDescription;
660  myDescription << Explicit3Stage3rdOrderHeun_name() << "\n"
661  << "Solving Ordinary Differential Equations I:\n"
662  << "Nonstiff Problems, 2nd Revised Edition\n"
663  << "E. Hairer, S.P. Norsett, G. Wanner\n"
664  << "Table 1.1, pg 135\n"
665  << "c = [ 0 1/3 2/3 ]'\n"
666  << "A = [ 0 ] \n"
667  << " [ 1/3 0 ]\n"
668  << " [ 0 2/3 0 ]\n"
669  << "b = [ 1/4 0 3/4 ]'" << std::endl;
670  typedef ScalarTraits<Scalar> ST;
671  Scalar one = ST::one();
672  Scalar zero = ST::zero();
673  Scalar onethird = one/(3*one);
674  Scalar twothirds = 2*one/(3*one);
675  Scalar onefourth = one/(4*one);
676  Scalar threefourths = 3*one/(4*one);
677 
678  int myNumStages = 3;
679  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
680  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
681  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
682 
683  // Fill myA:
684  myA(0,0) = zero;
685  myA(0,1) = zero;
686  myA(0,2) = zero;
687 
688  myA(1,0) = onethird;
689  myA(1,1) = zero;
690  myA(1,2) = zero;
691 
692  myA(2,0) = zero;
693  myA(2,1) = twothirds;
694  myA(2,2) = zero;
695 
696  // Fill myb:
697  myb(0) = onefourth;
698  myb(1) = zero;
699  myb(2) = threefourths;
700 
701  // fill b_c_
702  myc(0) = zero;
703  myc(1) = onethird;
704  myc(2) = twothirds;
705 
706  this->setMyDescription(myDescription.str());
707  this->setMy_A(myA);
708  this->setMy_b(myb);
709  this->setMy_c(myc);
710  this->setMy_order(3);
711  }
712 };
713 
714 
715 template<class Scalar>
716 class Explicit2Stage2ndOrderRunge_RKBT :
717  virtual public RKButcherTableauDefaultBase<Scalar>
718 {
719  public:
720  Explicit2Stage2ndOrderRunge_RKBT()
721  {
722  std::ostringstream myDescription;
723  myDescription << Explicit2Stage2ndOrderRunge_name() << "\n"
724  << "Also known as Explicit Midpoint\n"
725  << "Solving Ordinary Differential Equations I:\n"
726  << "Nonstiff Problems, 2nd Revised Edition\n"
727  << "E. Hairer, S.P. Norsett, G. Wanner\n"
728  << "Table 1.1, pg 135\n"
729  << "c = [ 0 1/2 ]'\n"
730  << "A = [ 0 ]\n"
731  << " [ 1/2 0 ]\n"
732  << "b = [ 0 1 ]'" << std::endl;
733  typedef ScalarTraits<Scalar> ST;
734  Scalar one = ST::one();
735  Scalar zero = ST::zero();
736  Scalar onehalf = ST::one()/(2*ST::one());
737 
738  int myNumStages = 2;
739  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
740  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
741  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
742 
743  // Fill myA:
744  myA(0,0) = zero;
745  myA(0,1) = zero;
746 
747  myA(1,0) = onehalf;
748  myA(1,1) = zero;
749 
750  // Fill myb:
751  myb(0) = zero;
752  myb(1) = one;
753 
754  // fill b_c_
755  myc(0) = zero;
756  myc(1) = onehalf;
757 
758  this->setMyDescription(myDescription.str());
759  this->setMy_A(myA);
760  this->setMy_b(myb);
761  this->setMy_c(myc);
762  this->setMy_order(2);
763  }
764 };
765 
766 
767 template<class Scalar>
768 class ExplicitTrapezoidal_RKBT :
769  virtual public RKButcherTableauDefaultBase<Scalar>
770 {
771  public:
772  ExplicitTrapezoidal_RKBT()
773  {
774  std::ostringstream myDescription;
775  myDescription << ExplicitTrapezoidal_name() << "\n"
776  << "c = [ 0 1 ]'\n"
777  << "A = [ 0 ]\n"
778  << " [ 1 0 ]\n"
779  << "b = [ 1/2 1/2 ]'" << std::endl;
780  typedef ScalarTraits<Scalar> ST;
781  Scalar one = ST::one();
782  Scalar zero = ST::zero();
783  Scalar onehalf = ST::one()/(2*ST::one());
784 
785  int myNumStages = 2;
786  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
787  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
788  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
789 
790  // Fill myA:
791  myA(0,0) = zero;
792  myA(0,1) = zero;
793 
794  myA(1,0) = one;
795  myA(1,1) = zero;
796 
797  // Fill myb:
798  myb(0) = onehalf;
799  myb(1) = onehalf;
800 
801  // fill b_c_
802  myc(0) = zero;
803  myc(1) = one;
804 
805  this->setMyDescription(myDescription.str());
806  this->setMy_A(myA);
807  this->setMy_b(myb);
808  this->setMy_c(myc);
809  this->setMy_order(2);
810  }
811 };
812 
813 
814 template<class Scalar>
815 class SDIRK2Stage2ndOrder_RKBT :
816  virtual public RKButcherTableauDefaultBase<Scalar>
817 {
818  public:
819  SDIRK2Stage2ndOrder_RKBT()
820  {
821  std::ostringstream myDescription;
822  myDescription << SDIRK2Stage2ndOrder_name() << "\n"
823  << "Computer Methods for ODEs and DAEs\n"
824  << "U. M. Ascher and L. R. Petzold\n"
825  << "p. 106\n"
826  << "gamma = (2+-sqrt(2))/2\n"
827  << "c = [ gamma 1 ]'\n"
828  << "A = [ gamma 0 ]\n"
829  << " [ 1-gamma gamma ]\n"
830  << "b = [ 1-gamma gamma ]'" << std::endl;
831 
832  this->setMyDescription(myDescription.str());
833  typedef ScalarTraits<Scalar> ST;
834  Scalar one = ST::one();
835  gamma_default_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
836  gamma_ = gamma_default_;
837  this->setupData();
838 
839  RCP<ParameterList> validPL = Teuchos::parameterList();
840  validPL->set("Description","",this->getMyDescription());
841  validPL->set<double>("gamma",gamma_default_,
842  "The default value is gamma = (2-sqrt(2))/2. "
843  "This will produce an L-stable 2nd order method with the stage "
844  "times within the timestep. Other values of gamma will still "
845  "produce an L-stable scheme, but will only be 1st order accurate.");
846  Teuchos::setupVerboseObjectSublist(&*validPL);
847  this->setMyValidParameterList(validPL);
848  }
849  void setupData()
850  {
851  typedef ScalarTraits<Scalar> ST;
852  int myNumStages = 2;
853  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
854  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
855  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
856  Scalar one = ST::one();
857  Scalar zero = ST::zero();
858  myA(0,0) = gamma_;
859  myA(0,1) = zero;
860  myA(1,0) = as<Scalar>( one - gamma_ );
861  myA(1,1) = gamma_;
862  myb(0) = as<Scalar>( one - gamma_ );
863  myb(1) = gamma_;
864  myc(0) = gamma_;
865  myc(1) = one;
866 
867  this->setMy_A(myA);
868  this->setMy_b(myb);
869  this->setMy_c(myc);
870  this->setMy_order(2);
871  }
872  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
873  {
874  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
875  paramList->validateParameters(*this->getValidParameters());
876  Teuchos::readVerboseObjectSublist(&*paramList,this);
877  gamma_ = paramList->get<double>("gamma",gamma_default_);
878  this->setupData();
879  this->setMyParamList(paramList);
880  }
881  private:
882  Scalar gamma_default_;
883  Scalar gamma_;
884 };
885 
886 
887 // 04/07/09 tscoffe: I verified manually that the Convergence Testing passes
888 // with gamma_default_ = -1.
889 template<class Scalar>
890 class SDIRK2Stage3rdOrder_RKBT :
891  virtual public RKButcherTableauDefaultBase<Scalar>
892 {
893  public:
894  SDIRK2Stage3rdOrder_RKBT()
895  {
896  std::ostringstream myDescription;
897  myDescription << SDIRK2Stage3rdOrder_name() << "\n"
898  << "Solving Ordinary Differential Equations I:\n"
899  << "Nonstiff Problems, 2nd Revised Edition\n"
900  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
901  << "Table 7.2, pg 207\n"
902  << "gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n"
903  << "gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n"
904  << "c = [ gamma 1-gamma ]'\n"
905  << "A = [ gamma 0 ]\n"
906  << " [ 1-2*gamma gamma ]\n"
907  << "b = [ 1/2 1/2 ]'" << std::endl;
908 
909  this->setMyDescription(myDescription.str());
910  thirdOrderAStable_default_ = true;
911  secondOrderLStable_default_ = false;
912  thirdOrderAStable_ = true;
913  secondOrderLStable_ = false;
914  typedef ScalarTraits<Scalar> ST;
915  Scalar one = ST::one();
916  gamma_default_ = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
917  gamma_ = gamma_default_;
918  this->setupData();
919 
920  RCP<ParameterList> validPL = Teuchos::parameterList();
921  validPL->set("Description","",this->getMyDescription());
922  validPL->set("3rd Order A-stable",thirdOrderAStable_default_,
923  "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain "
924  "a 3rd order A-stable scheme. '3rd Order A-stable' and "
925  "'2nd Order L-stable' can not both be true.");
926  validPL->set("2nd Order L-stable",secondOrderLStable_default_,
927  "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain "
928  "a 2nd order L-stable scheme. '3rd Order A-stable' and "
929  "'2nd Order L-stable' can not both be true.");
930  validPL->set("gamma",gamma_default_,
931  "If both '3rd Order A-stable' and '2nd Order L-stable' "
932  "are false, gamma will be used. The default value is the "
933  "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
934 
935  Teuchos::setupVerboseObjectSublist(&*validPL);
936  this->setMyValidParameterList(validPL);
937  }
938  void setupData()
939  {
940  typedef ScalarTraits<Scalar> ST;
941  int myNumStages = 2;
942  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
943  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
944  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
945  Scalar one = ST::one();
946  Scalar zero = ST::zero();
947  Scalar gamma;
948  if (thirdOrderAStable_)
949  gamma = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
950  else if (secondOrderLStable_)
951  gamma = as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
952  else
953  gamma = gamma_;
954  myA(0,0) = gamma;
955  myA(0,1) = zero;
956  myA(1,0) = as<Scalar>( one - 2*gamma );
957  myA(1,1) = gamma;
958  myb(0) = as<Scalar>( one/(2*one) );
959  myb(1) = as<Scalar>( one/(2*one) );
960  myc(0) = gamma;
961  myc(1) = as<Scalar>( one - gamma );
962 
963  this->setMy_A(myA);
964  this->setMy_b(myb);
965  this->setMy_c(myc);
966  this->setMy_order(3);
967  }
968  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
969  {
970  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
971  paramList->validateParameters(*this->getValidParameters());
972  Teuchos::readVerboseObjectSublist(&*paramList,this);
973  thirdOrderAStable_ = paramList->get("3rd Order A-stable",
974  thirdOrderAStable_default_);
975  secondOrderLStable_ = paramList->get("2nd Order L-stable",
976  secondOrderLStable_default_);
977  TEUCHOS_TEST_FOR_EXCEPTION(
978  thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
979  "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
980  gamma_ = paramList->get("gamma",gamma_default_);
981 
982  this->setupData();
983  this->setMyParamList(paramList);
984  }
985  private:
986  bool thirdOrderAStable_default_;
987  bool thirdOrderAStable_;
988  bool secondOrderLStable_default_;
989  bool secondOrderLStable_;
990  Scalar gamma_default_;
991  Scalar gamma_;
992 };
993 
994 
995 template<class Scalar>
996 class DIRK2Stage3rdOrder_RKBT :
997  virtual public RKButcherTableauDefaultBase<Scalar>
998 {
999  public:
1000  DIRK2Stage3rdOrder_RKBT()
1001  {
1002 
1003  std::ostringstream myDescription;
1004  myDescription << DIRK2Stage3rdOrder_name() << "\n"
1005  << "Hammer & Hollingsworth method\n"
1006  << "Solving Ordinary Differential Equations I:\n"
1007  << "Nonstiff Problems, 2nd Revised Edition\n"
1008  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1009  << "Table 7.1, pg 205\n"
1010  << "c = [ 0 2/3 ]'\n"
1011  << "A = [ 0 0 ]\n"
1012  << " [ 1/3 1/3 ]\n"
1013  << "b = [ 1/4 3/4 ]'" << std::endl;
1014  typedef ScalarTraits<Scalar> ST;
1015  int myNumStages = 2;
1016  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1017  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1018  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1019  Scalar one = ST::one();
1020  Scalar zero = ST::zero();
1021  myA(0,0) = zero;
1022  myA(0,1) = zero;
1023  myA(1,0) = as<Scalar>( one/(3*one) );
1024  myA(1,1) = as<Scalar>( one/(3*one) );
1025  myb(0) = as<Scalar>( one/(4*one) );
1026  myb(1) = as<Scalar>( 3*one/(4*one) );
1027  myc(0) = zero;
1028  myc(1) = as<Scalar>( 2*one/(3*one) );
1029  this->setMyDescription(myDescription.str());
1030  this->setMy_A(myA);
1031  this->setMy_b(myb);
1032  this->setMy_c(myc);
1033  this->setMy_order(3);
1034  }
1035 };
1036 
1037 
1038 template<class Scalar>
1039 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
1040  virtual public RKButcherTableauDefaultBase<Scalar>
1041 {
1042  public:
1043  Implicit3Stage6thOrderKuntzmannButcher_RKBT()
1044  {
1045 
1046  std::ostringstream myDescription;
1047  myDescription << Implicit3Stage6thOrderKuntzmannButcher_name() << "\n"
1048  << "Kuntzmann & Butcher method\n"
1049  << "Solving Ordinary Differential Equations I:\n"
1050  << "Nonstiff Problems, 2nd Revised Edition\n"
1051  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1052  << "Table 7.4, pg 209\n"
1053  << "c = [ 1/2-sqrt(15)/10 1/2 1/2-sqrt(15)/10 ]'\n"
1054  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1055  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1056  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1057  << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1058  typedef ScalarTraits<Scalar> ST;
1059  int myNumStages = 3;
1060  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1061  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1062  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1063  Scalar one = ST::one();
1064  myA(0,0) = as<Scalar>( 5*one/(36*one) );
1065  myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
1066  myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
1067  myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
1068  myA(1,1) = as<Scalar>( 2*one/(9*one) );
1069  myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
1070  myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
1071  myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
1072  myA(2,2) = as<Scalar>( 5*one/(36*one) );
1073  myb(0) = as<Scalar>( 5*one/(18*one) );
1074  myb(1) = as<Scalar>( 4*one/(9*one) );
1075  myb(2) = as<Scalar>( 5*one/(18*one) );
1076  myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
1077  myc(1) = as<Scalar>( one/(2*one) );
1078  myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
1079  this->setMyDescription(myDescription.str());
1080  this->setMy_A(myA);
1081  this->setMy_b(myb);
1082  this->setMy_c(myc);
1083  this->setMy_order(6);
1084  }
1085 };
1086 
1087 
1088 template<class Scalar>
1089 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
1090  virtual public RKButcherTableauDefaultBase<Scalar>
1091 {
1092  public:
1093  Implicit4Stage8thOrderKuntzmannButcher_RKBT()
1094  {
1095 
1096  std::ostringstream myDescription;
1097  myDescription << Implicit4Stage8thOrderKuntzmannButcher_name() << "\n"
1098  << "Kuntzmann & Butcher method\n"
1099  << "Solving Ordinary Differential Equations I:\n"
1100  << "Nonstiff Problems, 2nd Revised Edition\n"
1101  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1102  << "Table 7.5, pg 209\n"
1103  << "c = [ 1/2-w2 1/2-w2p 1/2+w2p 1/2+w2 ]'\n"
1104  << "A = [ w1 w1p-w3+w4p w1p-w3-w4p w1-w5 ]\n"
1105  << " [ w1-w3p+w4 w1p w1p-w5p w1-w3p-w4 ]\n"
1106  << " [ w1+w3p+w4 w1p+w5p w1p w1+w3p-w4 ]\n"
1107  << " [ w1+w5 w1p+w3+w4p w1p+w3-w4p w1 ]\n"
1108  << "b = [ 2*w1 2*w1p 2*w1p 2*w1 ]'\n"
1109  << "w1 = 1/8-sqrt(30)/144\n"
1110  << "w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n"
1111  << "w3 = w2*(1/6+sqrt(30)/24)\n"
1112  << "w4 = w2*(1/21+5*sqrt(30)/168)\n"
1113  << "w5 = w2-2*w3\n"
1114  << "w1p = 1/8+sqrt(30)/144\n"
1115  << "w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n"
1116  << "w3p = w2*(1/6-sqrt(30)/24)\n"
1117  << "w4p = w2*(1/21-5*sqrt(30)/168)\n"
1118  << "w5p = w2p-2*w3p" << std::endl;
1119  typedef ScalarTraits<Scalar> ST;
1120  int myNumStages = 4;
1121  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1122  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1123  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1124  Scalar one = ST::one();
1125  Scalar onehalf = as<Scalar>( one/(2*one) );
1126  Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
1127  Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
1128  Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
1129  Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
1130  Scalar w5 = as<Scalar>( w2-2*w3 );
1131  Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
1132  Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
1133  Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
1134  Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
1135  Scalar w5p = as<Scalar>( w2p-2*w3p );
1136  myA(0,0) = w1;
1137  myA(0,1) = w1p-w3+w4p;
1138  myA(0,2) = w1p-w3-w4p;
1139  myA(0,3) = w1-w5;
1140  myA(1,0) = w1-w3p+w4;
1141  myA(1,1) = w1p;
1142  myA(1,2) = w1p-w5p;
1143  myA(1,3) = w1-w3p-w4;
1144  myA(2,0) = w1+w3p+w4;
1145  myA(2,1) = w1p+w5p;
1146  myA(2,2) = w1p;
1147  myA(2,3) = w1+w3p-w4;
1148  myA(3,0) = w1+w5;
1149  myA(3,1) = w1p+w3+w4p;
1150  myA(3,2) = w1p+w3-w4p;
1151  myA(3,3) = w1;
1152  myb(0) = 2*w1;
1153  myb(1) = 2*w1p;
1154  myb(2) = 2*w1p;
1155  myb(3) = 2*w1;
1156  myc(0) = onehalf - w2;
1157  myc(1) = onehalf - w2p;
1158  myc(2) = onehalf + w2p;
1159  myc(3) = onehalf + w2;
1160  this->setMyDescription(myDescription.str());
1161  this->setMy_A(myA);
1162  this->setMy_b(myb);
1163  this->setMy_c(myc);
1164  this->setMy_order(8);
1165  }
1166 };
1167 
1168 
1169 template<class Scalar>
1170 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
1171  virtual public RKButcherTableauDefaultBase<Scalar>
1172 {
1173  public:
1174  Implicit2Stage4thOrderHammerHollingsworth_RKBT()
1175  {
1176 
1177  std::ostringstream myDescription;
1178  myDescription << Implicit2Stage4thOrderHammerHollingsworth_name() << "\n"
1179  << "Hammer & Hollingsworth method\n"
1180  << "Solving Ordinary Differential Equations I:\n"
1181  << "Nonstiff Problems, 2nd Revised Edition\n"
1182  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1183  << "Table 7.3, pg 207\n"
1184  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1185  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1186  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1187  << "b = [ 1/2 1/2 ]'" << std::endl;
1188  typedef ScalarTraits<Scalar> ST;
1189  int myNumStages = 2;
1190  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1191  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1192  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1193  Scalar one = ST::one();
1194  Scalar onequarter = as<Scalar>( one/(4*one) );
1195  Scalar onehalf = as<Scalar>( one/(2*one) );
1196  myA(0,0) = onequarter;
1197  myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
1198  myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
1199  myA(1,1) = onequarter;
1200  myb(0) = onehalf;
1201  myb(1) = onehalf;
1202  myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
1203  myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
1204  this->setMyDescription(myDescription.str());
1205  this->setMy_A(myA);
1206  this->setMy_b(myb);
1207  this->setMy_c(myc);
1208  this->setMy_order(4);
1209  }
1210 };
1211 
1212 
1213 template<class Scalar>
1214 class IRK1StageTheta_RKBT :
1215  virtual public RKButcherTableauDefaultBase<Scalar>
1216 {
1217  public:
1218  IRK1StageTheta_RKBT()
1219  {
1220 
1221  std::ostringstream myDescription;
1222  myDescription << IRK1StageTheta_name() << "\n"
1223  << "Non-standard finite-difference methods\n"
1224  << "in dynamical systems, P. Kama,\n"
1225  << "Dissertation, University of Pretoria, pg. 49.\n"
1226  << "Comment: Generalized Implicit Midpoint Method\n"
1227  << "c = [ theta ]'\n"
1228  << "A = [ theta ]\n"
1229  << "b = [ 1 ]'\n"
1230  << std::endl;
1231 
1232  this->setMyDescription(myDescription.str());
1233  typedef ScalarTraits<Scalar> ST;
1234  theta_default_ = ST::one()/(2*ST::one());
1235  theta_ = theta_default_;
1236  this->setupData();
1237 
1238  RCP<ParameterList> validPL = Teuchos::parameterList();
1239  validPL->set("Description","",this->getMyDescription());
1240  validPL->set<double>("theta",theta_default_,
1241  "Valid values are 0 <= theta <= 1, where theta = 0 "
1242  "implies Forward Euler, theta = 1/2 implies midpoint "
1243  "method, and theta = 1 implies Backward Euler. "
1244  "For theta != 1/2, this method is first-order accurate, "
1245  "and with theta = 1/2, it is second-order accurate. "
1246  "This method is A-stable, but becomes L-stable with theta=1.");
1247  Teuchos::setupVerboseObjectSublist(&*validPL);
1248  this->setMyValidParameterList(validPL);
1249  }
1250 
1251  void setupData()
1252  {
1253  typedef ScalarTraits<Scalar> ST;
1254  int myNumStages = 1;
1255  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1256  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1257  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1258  myA(0,0) = theta_;
1259  myb(0) = ST::one();
1260  myc(0) = theta_;
1261  this->setMy_A(myA);
1262  this->setMy_b(myb);
1263  this->setMy_c(myc);
1264  if (theta_ == theta_default_) this->setMy_order(2);
1265  else this->setMy_order(1);
1266  }
1267 
1268  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1269  {
1270  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1271  paramList->validateParameters(*this->getValidParameters());
1272  Teuchos::readVerboseObjectSublist(&*paramList,this);
1273  theta_ = paramList->get<double>("theta",theta_default_);
1274  this->setupData();
1275  this->setMyParamList(paramList);
1276  }
1277  private:
1278  Scalar theta_default_;
1279  Scalar theta_;
1280 };
1281 
1282 
1283 template<class Scalar>
1284 class IRK2StageTheta_RKBT :
1285  virtual public RKButcherTableauDefaultBase<Scalar>
1286 {
1287  public:
1288  IRK2StageTheta_RKBT()
1289  {
1290  std::ostringstream myDescription;
1291  myDescription << IRK2StageTheta_name() << "\n"
1292  << "Computer Methods for ODEs and DAEs\n"
1293  << "U. M. Ascher and L. R. Petzold\n"
1294  << "p. 113\n"
1295  << "c = [ 0 1 ]'\n"
1296  << "A = [ 0 0 ]\n"
1297  << " [ 1-theta theta ]\n"
1298  << "b = [ 1-theta theta ]'\n"
1299  << std::endl;
1300 
1301  this->setMyDescription(myDescription.str());
1302  typedef ScalarTraits<Scalar> ST;
1303  theta_default_ = ST::one()/(2*ST::one());
1304  theta_ = theta_default_;
1305  this->setupData();
1306 
1307  RCP<ParameterList> validPL = Teuchos::parameterList();
1308  validPL->set("Description","",this->getMyDescription());
1309  validPL->set<double>("theta",theta_default_,
1310  "Valid values are 0 < theta <= 1, where theta = 0 "
1311  "implies Forward Euler, theta = 1/2 implies trapezoidal "
1312  "method, and theta = 1 implies Backward Euler. "
1313  "For theta != 1/2, this method is first-order accurate, "
1314  "and with theta = 1/2, it is second-order accurate. "
1315  "This method is A-stable, but becomes L-stable with theta=1.");
1316  Teuchos::setupVerboseObjectSublist(&*validPL);
1317  this->setMyValidParameterList(validPL);
1318  }
1319  void setupData()
1320  {
1321  typedef ScalarTraits<Scalar> ST;
1322  int myNumStages = 2;
1323  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1324  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1325  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1326  Scalar one = ST::one();
1327  Scalar zero = ST::zero();
1328  myA(0,0) = zero;
1329  myA(0,1) = zero;
1330  myA(1,0) = as<Scalar>( one - theta_ );
1331  myA(1,1) = theta_;
1332  myb(0) = as<Scalar>( one - theta_ );
1333  myb(1) = theta_;
1334  myc(0) = theta_;
1335  myc(1) = one;
1336 
1337  this->setMy_A(myA);
1338  this->setMy_b(myb);
1339  this->setMy_c(myc);
1340  TEUCHOS_TEST_FOR_EXCEPTION(
1341  theta_ == zero, std::logic_error,
1342  "'theta' can not be zero, as it makes this IRK stepper explicit.");
1343  if (theta_ == theta_default_) this->setMy_order(2);
1344  else this->setMy_order(1);
1345  }
1346  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1347  {
1348  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1349  paramList->validateParameters(*this->getValidParameters());
1350  Teuchos::readVerboseObjectSublist(&*paramList,this);
1351  theta_ = paramList->get<double>("theta",theta_default_);
1352  this->setupData();
1353  this->setMyParamList(paramList);
1354  }
1355  private:
1356  Scalar theta_default_;
1357  Scalar theta_;
1358 };
1359 
1360 
1361 template<class Scalar>
1362 class Implicit1Stage2ndOrderGauss_RKBT :
1363  virtual public RKButcherTableauDefaultBase<Scalar>
1364 {
1365  public:
1366  Implicit1Stage2ndOrderGauss_RKBT()
1367  {
1368 
1369  std::ostringstream myDescription;
1370  myDescription << Implicit1Stage2ndOrderGauss_name() << "\n"
1371  << "A-stable\n"
1372  << "Solving Ordinary Differential Equations II:\n"
1373  << "Stiff and Differential-Algebraic Problems,\n"
1374  << "2nd Revised Edition\n"
1375  << "E. Hairer and G. Wanner\n"
1376  << "Table 5.2, pg 72\n"
1377  << "Also: Implicit midpoint rule\n"
1378  << "Solving Ordinary Differential Equations I:\n"
1379  << "Nonstiff Problems, 2nd Revised Edition\n"
1380  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1381  << "Table 7.1, pg 205\n"
1382  << "c = [ 1/2 ]'\n"
1383  << "A = [ 1/2 ]\n"
1384  << "b = [ 1 ]'" << std::endl;
1385  typedef ScalarTraits<Scalar> ST;
1386  int myNumStages = 1;
1387  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1388  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1389  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1390  Scalar onehalf = ST::one()/(2*ST::one());
1391  Scalar one = ST::one();
1392  myA(0,0) = onehalf;
1393  myb(0) = one;
1394  myc(0) = onehalf;
1395  this->setMyDescription(myDescription.str());
1396  this->setMy_A(myA);
1397  this->setMy_b(myb);
1398  this->setMy_c(myc);
1399  this->setMy_order(2);
1400  }
1401 };
1402 
1403 
1404 template<class Scalar>
1405 class Implicit2Stage4thOrderGauss_RKBT :
1406  virtual public RKButcherTableauDefaultBase<Scalar>
1407 {
1408  public:
1409  Implicit2Stage4thOrderGauss_RKBT()
1410  {
1411 
1412  std::ostringstream myDescription;
1413  myDescription << Implicit2Stage4thOrderGauss_name() << "\n"
1414  << "A-stable\n"
1415  << "Solving Ordinary Differential Equations II:\n"
1416  << "Stiff and Differential-Algebraic Problems,\n"
1417  << "2nd Revised Edition\n"
1418  << "E. Hairer and G. Wanner\n"
1419  << "Table 5.2, pg 72\n"
1420  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1421  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1422  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1423  << "b = [ 1/2 1/2 ]'" << std::endl;
1424  typedef ScalarTraits<Scalar> ST;
1425  int myNumStages = 2;
1426  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1427  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1428  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1429  Scalar one = ST::one();
1430  Scalar onehalf = as<Scalar>(one/(2*one));
1431  Scalar three = as<Scalar>(3*one);
1432  Scalar six = as<Scalar>(6*one);
1433  Scalar onefourth = as<Scalar>(one/(4*one));
1434  Scalar alpha = ST::squareroot(three)/six;
1435 
1436  myA(0,0) = onefourth;
1437  myA(0,1) = onefourth-alpha;
1438  myA(1,0) = onefourth+alpha;
1439  myA(1,1) = onefourth;
1440  myb(0) = onehalf;
1441  myb(1) = onehalf;
1442  myc(0) = onehalf-alpha;
1443  myc(1) = onehalf+alpha;
1444  this->setMyDescription(myDescription.str());
1445  this->setMy_A(myA);
1446  this->setMy_b(myb);
1447  this->setMy_c(myc);
1448  this->setMy_order(4);
1449  }
1450 };
1451 
1452 
1453 template<class Scalar>
1454 class Implicit3Stage6thOrderGauss_RKBT :
1455  virtual public RKButcherTableauDefaultBase<Scalar>
1456 {
1457  public:
1458  Implicit3Stage6thOrderGauss_RKBT()
1459  {
1460 
1461  std::ostringstream myDescription;
1462  myDescription << Implicit3Stage6thOrderGauss_name() << "\n"
1463  << "A-stable\n"
1464  << "Solving Ordinary Differential Equations II:\n"
1465  << "Stiff and Differential-Algebraic Problems,\n"
1466  << "2nd Revised Edition\n"
1467  << "E. Hairer and G. Wanner\n"
1468  << "Table 5.2, pg 72\n"
1469  << "c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n"
1470  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1471  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1472  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1473  << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1474  typedef ScalarTraits<Scalar> ST;
1475  int myNumStages = 3;
1476  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1477  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1478  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1479  Scalar one = ST::one();
1480  Scalar ten = as<Scalar>(10*one);
1481  Scalar fifteen = as<Scalar>(15*one);
1482  Scalar twentyfour = as<Scalar>(24*one);
1483  Scalar thirty = as<Scalar>(30*one);
1484  Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
1485  Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
1486  Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
1487  Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
1488 
1489  myA(0,0) = as<Scalar>(5*one/(36*one));
1490  myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
1491  myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
1492  myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
1493  myA(1,1) = as<Scalar>(2*one/(9*one));
1494  myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
1495  myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
1496  myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
1497  myA(2,2) = as<Scalar>(5*one/(36*one));
1498  myb(0) = as<Scalar>(5*one/(18*one));
1499  myb(1) = as<Scalar>(4*one/(9*one));
1500  myb(2) = as<Scalar>(5*one/(18*one));
1501  myc(0) = as<Scalar>(one/(2*one))-sqrt15over10;
1502  myc(1) = as<Scalar>(one/(2*one));
1503  myc(2) = as<Scalar>(one/(2*one))+sqrt15over10;
1504  this->setMyDescription(myDescription.str());
1505  this->setMy_A(myA);
1506  this->setMy_b(myb);
1507  this->setMy_c(myc);
1508  this->setMy_order(6);
1509  }
1510 };
1511 
1512 
1513 template<class Scalar>
1514 class Implicit1Stage1stOrderRadauA_RKBT :
1515  virtual public RKButcherTableauDefaultBase<Scalar>
1516 {
1517  public:
1518  Implicit1Stage1stOrderRadauA_RKBT()
1519  {
1520 
1521  std::ostringstream myDescription;
1522  myDescription << Implicit1Stage1stOrderRadauA_name() << "\n"
1523  << "A-stable\n"
1524  << "Solving Ordinary Differential Equations II:\n"
1525  << "Stiff and Differential-Algebraic Problems,\n"
1526  << "2nd Revised Edition\n"
1527  << "E. Hairer and G. Wanner\n"
1528  << "Table 5.3, pg 73\n"
1529  << "c = [ 0 ]'\n"
1530  << "A = [ 1 ]\n"
1531  << "b = [ 1 ]'" << std::endl;
1532  typedef ScalarTraits<Scalar> ST;
1533  int myNumStages = 1;
1534  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1535  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1536  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1537  Scalar one = ST::one();
1538  Scalar zero = ST::zero();
1539  myA(0,0) = one;
1540  myb(0) = one;
1541  myc(0) = zero;
1542  this->setMyDescription(myDescription.str());
1543  this->setMy_A(myA);
1544  this->setMy_b(myb);
1545  this->setMy_c(myc);
1546  this->setMy_order(1);
1547  }
1548 };
1549 
1550 
1551 template<class Scalar>
1552 class Implicit2Stage3rdOrderRadauA_RKBT :
1553  virtual public RKButcherTableauDefaultBase<Scalar>
1554 {
1555  public:
1556  Implicit2Stage3rdOrderRadauA_RKBT()
1557  {
1558 
1559  std::ostringstream myDescription;
1560  myDescription << Implicit2Stage3rdOrderRadauA_name() << "\n"
1561  << "A-stable\n"
1562  << "Solving Ordinary Differential Equations II:\n"
1563  << "Stiff and Differential-Algebraic Problems,\n"
1564  << "2nd Revised Edition\n"
1565  << "E. Hairer and G. Wanner\n"
1566  << "Table 5.3, pg 73\n"
1567  << "c = [ 0 2/3 ]'\n"
1568  << "A = [ 1/4 -1/4 ]\n"
1569  << " [ 1/4 5/12 ]\n"
1570  << "b = [ 1/4 3/4 ]'" << std::endl;
1571  typedef ScalarTraits<Scalar> ST;
1572  int myNumStages = 2;
1573  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1574  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1575  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1576  Scalar zero = ST::zero();
1577  Scalar one = ST::one();
1578  myA(0,0) = as<Scalar>(one/(4*one));
1579  myA(0,1) = as<Scalar>(-one/(4*one));
1580  myA(1,0) = as<Scalar>(one/(4*one));
1581  myA(1,1) = as<Scalar>(5*one/(12*one));
1582  myb(0) = as<Scalar>(one/(4*one));
1583  myb(1) = as<Scalar>(3*one/(4*one));
1584  myc(0) = zero;
1585  myc(1) = as<Scalar>(2*one/(3*one));
1586  this->setMyDescription(myDescription.str());
1587  this->setMy_A(myA);
1588  this->setMy_b(myb);
1589  this->setMy_c(myc);
1590  this->setMy_order(3);
1591  }
1592 };
1593 
1594 
1595 template<class Scalar>
1596 class Implicit3Stage5thOrderRadauA_RKBT :
1597  virtual public RKButcherTableauDefaultBase<Scalar>
1598 {
1599  public:
1600  Implicit3Stage5thOrderRadauA_RKBT()
1601  {
1602 
1603  std::ostringstream myDescription;
1604  myDescription << Implicit3Stage5thOrderRadauA_name() << "\n"
1605  << "A-stable\n"
1606  << "Solving Ordinary Differential Equations II:\n"
1607  << "Stiff and Differential-Algebraic Problems,\n"
1608  << "2nd Revised Edition\n"
1609  << "E. Hairer and G. Wanner\n"
1610  << "Table 5.4, pg 73\n"
1611  << "c = [ 0 (6-sqrt(6))/10 (6+sqrt(6))/10 ]'\n"
1612  << "A = [ 1/9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 ]\n"
1613  << " [ 1/9 (88+7*sqrt(6))/360 (88-43*sqrt(6))/360 ]\n"
1614  << " [ 1/9 (88+43*sqrt(6))/360 (88-7*sqrt(6))/360 ]\n"
1615  << "b = [ 1/9 (16+sqrt(6))/36 (16-sqrt(6))/36 ]'" << std::endl;
1616  typedef ScalarTraits<Scalar> ST;
1617  int myNumStages = 3;
1618  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1619  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1620  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1621  Scalar zero = ST::zero();
1622  Scalar one = ST::one();
1623  myA(0,0) = as<Scalar>(one/(9*one));
1624  myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
1625  myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
1626  myA(1,0) = as<Scalar>(one/(9*one));
1627  myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1628  myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
1629  myA(2,0) = as<Scalar>(one/(9*one));
1630  myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
1631  myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1632  myb(0) = as<Scalar>(one/(9*one));
1633  myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1634  myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1635  myc(0) = zero;
1636  myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
1637  myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
1638  this->setMyDescription(myDescription.str());
1639  this->setMy_A(myA);
1640  this->setMy_b(myb);
1641  this->setMy_c(myc);
1642  this->setMy_order(5);
1643  }
1644 };
1645 
1646 
1647 template<class Scalar>
1648 class Implicit1Stage1stOrderRadauB_RKBT :
1649  virtual public RKButcherTableauDefaultBase<Scalar>
1650 {
1651  public:
1652  Implicit1Stage1stOrderRadauB_RKBT()
1653  {
1654 
1655  std::ostringstream myDescription;
1656  myDescription << Implicit1Stage1stOrderRadauB_name() << "\n"
1657  << "A-stable\n"
1658  << "Solving Ordinary Differential Equations II:\n"
1659  << "Stiff and Differential-Algebraic Problems,\n"
1660  << "2nd Revised Edition\n"
1661  << "E. Hairer and G. Wanner\n"
1662  << "Table 5.5, pg 74\n"
1663  << "c = [ 1 ]'\n"
1664  << "A = [ 1 ]\n"
1665  << "b = [ 1 ]'" << std::endl;
1666  typedef ScalarTraits<Scalar> ST;
1667  int myNumStages = 1;
1668  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1669  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1670  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1671  Scalar one = ST::one();
1672  myA(0,0) = one;
1673  myb(0) = one;
1674  myc(0) = one;
1675  this->setMyDescription(myDescription.str());
1676  this->setMy_A(myA);
1677  this->setMy_b(myb);
1678  this->setMy_c(myc);
1679  this->setMy_order(1);
1680  }
1681 };
1682 
1683 
1684 template<class Scalar>
1685 class Implicit2Stage3rdOrderRadauB_RKBT :
1686  virtual public RKButcherTableauDefaultBase<Scalar>
1687 {
1688  public:
1689  Implicit2Stage3rdOrderRadauB_RKBT()
1690  {
1691 
1692  std::ostringstream myDescription;
1693  myDescription << Implicit2Stage3rdOrderRadauB_name() << "\n"
1694  << "A-stable\n"
1695  << "Solving Ordinary Differential Equations II:\n"
1696  << "Stiff and Differential-Algebraic Problems,\n"
1697  << "2nd Revised Edition\n"
1698  << "E. Hairer and G. Wanner\n"
1699  << "Table 5.5, pg 74\n"
1700  << "c = [ 1/3 1 ]'\n"
1701  << "A = [ 5/12 -1/12 ]\n"
1702  << " [ 3/4 1/4 ]\n"
1703  << "b = [ 3/4 1/4 ]'" << std::endl;
1704  typedef ScalarTraits<Scalar> ST;
1705  int myNumStages = 2;
1706  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1707  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1708  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1709  Scalar one = ST::one();
1710  myA(0,0) = as<Scalar>( 5*one/(12*one) );
1711  myA(0,1) = as<Scalar>( -one/(12*one) );
1712  myA(1,0) = as<Scalar>( 3*one/(4*one) );
1713  myA(1,1) = as<Scalar>( one/(4*one) );
1714  myb(0) = as<Scalar>( 3*one/(4*one) );
1715  myb(1) = as<Scalar>( one/(4*one) );
1716  myc(0) = as<Scalar>( one/(3*one) );
1717  myc(1) = one;
1718  this->setMyDescription(myDescription.str());
1719  this->setMy_A(myA);
1720  this->setMy_b(myb);
1721  this->setMy_c(myc);
1722  this->setMy_order(3);
1723  }
1724 };
1725 
1726 
1727 template<class Scalar>
1728 class Implicit3Stage5thOrderRadauB_RKBT :
1729  virtual public RKButcherTableauDefaultBase<Scalar>
1730 {
1731  public:
1732  Implicit3Stage5thOrderRadauB_RKBT()
1733  {
1734 
1735  std::ostringstream myDescription;
1736  myDescription << Implicit3Stage5thOrderRadauB_name() << "\n"
1737  << "A-stable\n"
1738  << "Solving Ordinary Differential Equations II:\n"
1739  << "Stiff and Differential-Algebraic Problems,\n"
1740  << "2nd Revised Edition\n"
1741  << "E. Hairer and G. Wanner\n"
1742  << "Table 5.6, pg 74\n"
1743  << "c = [ (4-sqrt(6))/10 (4+sqrt(6))/10 1 ]'\n"
1744  << "A = [ A1 A2 A3 ]\n"
1745  << " A1 = [ (88-7*sqrt(6))/360 ]\n"
1746  << " [ (296+169*sqrt(6))/1800 ]\n"
1747  << " [ (16-sqrt(6))/36 ]\n"
1748  << " A2 = [ (296-169*sqrt(6))/1800 ]\n"
1749  << " [ (88+7*sqrt(6))/360 ]\n"
1750  << " [ (16+sqrt(6))/36 ]\n"
1751  << " A3 = [ (-2+3*sqrt(6))/225 ]\n"
1752  << " [ (-2-3*sqrt(6))/225 ]\n"
1753  << " [ 1/9 ]\n"
1754  << "b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'"
1755  << std::endl;
1756  typedef ScalarTraits<Scalar> ST;
1757  int myNumStages = 3;
1758  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1759  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1760  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1761  Scalar one = ST::one();
1762  myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1763  myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
1764  myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
1765  myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
1766  myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1767  myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
1768  myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1769  myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1770  myA(2,2) = as<Scalar>( one/(9*one) );
1771  myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1772  myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1773  myb(2) = as<Scalar>( one/(9*one) );
1774  myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
1775  myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
1776  myc(2) = one;
1777  this->setMyDescription(myDescription.str());
1778  this->setMy_A(myA);
1779  this->setMy_b(myb);
1780  this->setMy_c(myc);
1781  this->setMy_order(5);
1782  }
1783 };
1784 
1785 
1786 template<class Scalar>
1787 class Implicit2Stage2ndOrderLobattoA_RKBT :
1788  virtual public RKButcherTableauDefaultBase<Scalar>
1789 {
1790  public:
1791  Implicit2Stage2ndOrderLobattoA_RKBT()
1792  {
1793 
1794  std::ostringstream myDescription;
1795  myDescription << Implicit2Stage2ndOrderLobattoA_name() << "\n"
1796  << "A-stable\n"
1797  << "Solving Ordinary Differential Equations II:\n"
1798  << "Stiff and Differential-Algebraic Problems,\n"
1799  << "2nd Revised Edition\n"
1800  << "E. Hairer and G. Wanner\n"
1801  << "Table 5.7, pg 75\n"
1802  << "c = [ 0 1 ]'\n"
1803  << "A = [ 0 0 ]\n"
1804  << " [ 1/2 1/2 ]\n"
1805  << "b = [ 1/2 1/2 ]'" << std::endl;
1806  typedef ScalarTraits<Scalar> ST;
1807  int myNumStages = 2;
1808  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1809  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1810  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1811  Scalar zero = ST::zero();
1812  Scalar one = ST::one();
1813  myA(0,0) = zero;
1814  myA(0,1) = zero;
1815  myA(1,0) = as<Scalar>( one/(2*one) );
1816  myA(1,1) = as<Scalar>( one/(2*one) );
1817  myb(0) = as<Scalar>( one/(2*one) );
1818  myb(1) = as<Scalar>( one/(2*one) );
1819  myc(0) = zero;
1820  myc(1) = one;
1821  this->setMyDescription(myDescription.str());
1822  this->setMy_A(myA);
1823  this->setMy_b(myb);
1824  this->setMy_c(myc);
1825  this->setMy_order(2);
1826  }
1827 };
1828 
1829 
1830 template<class Scalar>
1831 class Implicit3Stage4thOrderLobattoA_RKBT :
1832  virtual public RKButcherTableauDefaultBase<Scalar>
1833 {
1834  public:
1835  Implicit3Stage4thOrderLobattoA_RKBT()
1836  {
1837 
1838  std::ostringstream myDescription;
1839  myDescription << Implicit3Stage4thOrderLobattoA_name() << "\n"
1840  << "A-stable\n"
1841  << "Solving Ordinary Differential Equations II:\n"
1842  << "Stiff and Differential-Algebraic Problems,\n"
1843  << "2nd Revised Edition\n"
1844  << "E. Hairer and G. Wanner\n"
1845  << "Table 5.7, pg 75\n"
1846  << "c = [ 0 1/2 1 ]'\n"
1847  << "A = [ 0 0 0 ]\n"
1848  << " [ 5/24 1/3 -1/24 ]\n"
1849  << " [ 1/6 2/3 1/6 ]\n"
1850  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
1851  typedef ScalarTraits<Scalar> ST;
1852  int myNumStages = 3;
1853  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1854  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1855  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1856  Scalar zero = ST::zero();
1857  Scalar one = ST::one();
1858  myA(0,0) = zero;
1859  myA(0,1) = zero;
1860  myA(0,2) = zero;
1861  myA(1,0) = as<Scalar>( (5*one)/(24*one) );
1862  myA(1,1) = as<Scalar>( (one)/(3*one) );
1863  myA(1,2) = as<Scalar>( (-one)/(24*one) );
1864  myA(2,0) = as<Scalar>( (one)/(6*one) );
1865  myA(2,1) = as<Scalar>( (2*one)/(3*one) );
1866  myA(2,2) = as<Scalar>( (1*one)/(6*one) );
1867  myb(0) = as<Scalar>( (one)/(6*one) );
1868  myb(1) = as<Scalar>( (2*one)/(3*one) );
1869  myb(2) = as<Scalar>( (1*one)/(6*one) );
1870  myc(0) = zero;
1871  myc(1) = as<Scalar>( one/(2*one) );
1872  myc(2) = one;
1873  this->setMyDescription(myDescription.str());
1874  this->setMy_A(myA);
1875  this->setMy_b(myb);
1876  this->setMy_c(myc);
1877  this->setMy_order(4);
1878  }
1879 };
1880 
1881 
1882 template<class Scalar>
1883 class Implicit4Stage6thOrderLobattoA_RKBT :
1884  virtual public RKButcherTableauDefaultBase<Scalar>
1885 {
1886  public:
1887  Implicit4Stage6thOrderLobattoA_RKBT()
1888  {
1889 
1890  using Teuchos::as;
1891  typedef Teuchos::ScalarTraits<Scalar> ST;
1892 
1893  std::ostringstream myDescription;
1894  myDescription << Implicit4Stage6thOrderLobattoA_name() << "\n"
1895  << "A-stable\n"
1896  << "Solving Ordinary Differential Equations II:\n"
1897  << "Stiff and Differential-Algebraic Problems,\n"
1898  << "2nd Revised Edition\n"
1899  << "E. Hairer and G. Wanner\n"
1900  << "Table 5.8, pg 75\n"
1901  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
1902  << "A = [ A1 A2 A3 A4 ]\n"
1903  << " A1 = [ 0 ]\n"
1904  << " [ (11+sqrt(5)/120 ]\n"
1905  << " [ (11-sqrt(5)/120 ]\n"
1906  << " [ 1/12 ]\n"
1907  << " A2 = [ 0 ]\n"
1908  << " [ (25-sqrt(5))/120 ]\n"
1909  << " [ (25+13*sqrt(5))/120 ]\n"
1910  << " [ 5/12 ]\n"
1911  << " A3 = [ 0 ]\n"
1912  << " [ (25-13*sqrt(5))/120 ]\n"
1913  << " [ (25+sqrt(5))/120 ]\n"
1914  << " [ 5/12 ]\n"
1915  << " A4 = [ 0 ]\n"
1916  << " [ (-1+sqrt(5))/120 ]\n"
1917  << " [ (-1-sqrt(5))/120 ]\n"
1918  << " [ 1/12 ]\n"
1919  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
1920  typedef ScalarTraits<Scalar> ST;
1921  int myNumStages = 4;
1922  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1923  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1924  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1925  Scalar zero = ST::zero();
1926  Scalar one = ST::one();
1927  myA(0,0) = zero;
1928  myA(0,1) = zero;
1929  myA(0,2) = zero;
1930  myA(0,3) = zero;
1931  myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
1932  myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
1933  myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
1934  myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
1935  myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
1936  myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
1937  myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
1938  myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
1939  myA(3,0) = as<Scalar>( one/(12*one) );
1940  myA(3,1) = as<Scalar>( 5*one/(12*one) );
1941  myA(3,2) = as<Scalar>( 5*one/(12*one) );
1942  myA(3,3) = as<Scalar>( one/(12*one) );
1943  myb(0) = as<Scalar>( one/(12*one) );
1944  myb(1) = as<Scalar>( 5*one/(12*one) );
1945  myb(2) = as<Scalar>( 5*one/(12*one) );
1946  myb(3) = as<Scalar>( one/(12*one) );
1947  myc(0) = zero;
1948  myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
1949  myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
1950  myc(3) = one;
1951  this->setMyDescription(myDescription.str());
1952  this->setMy_A(myA);
1953  this->setMy_b(myb);
1954  this->setMy_c(myc);
1955  this->setMy_order(6);
1956  }
1957 };
1958 
1959 
1960 template<class Scalar>
1961 class Implicit2Stage2ndOrderLobattoB_RKBT :
1962  virtual public RKButcherTableauDefaultBase<Scalar>
1963 {
1964  public:
1965  Implicit2Stage2ndOrderLobattoB_RKBT()
1966  {
1967 
1968  std::ostringstream myDescription;
1969  myDescription << Implicit2Stage2ndOrderLobattoB_name() << "\n"
1970  << "A-stable\n"
1971  << "Solving Ordinary Differential Equations II:\n"
1972  << "Stiff and Differential-Algebraic Problems,\n"
1973  << "2nd Revised Edition\n"
1974  << "E. Hairer and G. Wanner\n"
1975  << "Table 5.9, pg 76\n"
1976  << "c = [ 0 1 ]'\n"
1977  << "A = [ 1/2 0 ]\n"
1978  << " [ 1/2 0 ]\n"
1979  << "b = [ 1/2 1/2 ]'" << std::endl;
1980  typedef ScalarTraits<Scalar> ST;
1981  int myNumStages = 2;
1982  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1983  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1984  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1985  Scalar zero = ST::zero();
1986  Scalar one = ST::one();
1987  myA(0,0) = as<Scalar>( one/(2*one) );
1988  myA(0,1) = zero;
1989  myA(1,0) = as<Scalar>( one/(2*one) );
1990  myA(1,1) = zero;
1991  myb(0) = as<Scalar>( one/(2*one) );
1992  myb(1) = as<Scalar>( one/(2*one) );
1993  myc(0) = zero;
1994  myc(1) = one;
1995  this->setMyDescription(myDescription.str());
1996  this->setMy_A(myA);
1997  this->setMy_b(myb);
1998  this->setMy_c(myc);
1999  this->setMy_order(2);
2000  }
2001 };
2002 
2003 
2004 template<class Scalar>
2005 class Implicit3Stage4thOrderLobattoB_RKBT :
2006  virtual public RKButcherTableauDefaultBase<Scalar>
2007 {
2008  public:
2009  Implicit3Stage4thOrderLobattoB_RKBT()
2010  {
2011 
2012  std::ostringstream myDescription;
2013  myDescription << Implicit3Stage4thOrderLobattoB_name() << "\n"
2014  << "A-stable\n"
2015  << "Solving Ordinary Differential Equations II:\n"
2016  << "Stiff and Differential-Algebraic Problems,\n"
2017  << "2nd Revised Edition\n"
2018  << "E. Hairer and G. Wanner\n"
2019  << "Table 5.9, pg 76\n"
2020  << "c = [ 0 1/2 1 ]'\n"
2021  << "A = [ 1/6 -1/6 0 ]\n"
2022  << " [ 1/6 1/3 0 ]\n"
2023  << " [ 1/6 5/6 0 ]\n"
2024  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2025  typedef ScalarTraits<Scalar> ST;
2026  int myNumStages = 3;
2027  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2028  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2029  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2030  Scalar zero = ST::zero();
2031  Scalar one = ST::one();
2032  myA(0,0) = as<Scalar>( one/(6*one) );
2033  myA(0,1) = as<Scalar>( -one/(6*one) );
2034  myA(0,2) = zero;
2035  myA(1,0) = as<Scalar>( one/(6*one) );
2036  myA(1,1) = as<Scalar>( one/(3*one) );
2037  myA(1,2) = zero;
2038  myA(2,0) = as<Scalar>( one/(6*one) );
2039  myA(2,1) = as<Scalar>( 5*one/(6*one) );
2040  myA(2,2) = zero;
2041  myb(0) = as<Scalar>( one/(6*one) );
2042  myb(1) = as<Scalar>( 2*one/(3*one) );
2043  myb(2) = as<Scalar>( one/(6*one) );
2044  myc(0) = zero;
2045  myc(1) = as<Scalar>( one/(2*one) );
2046  myc(2) = one;
2047  this->setMyDescription(myDescription.str());
2048  this->setMy_A(myA);
2049  this->setMy_b(myb);
2050  this->setMy_c(myc);
2051  this->setMy_order(4);
2052  }
2053 };
2054 
2055 
2056 template<class Scalar>
2057 class Implicit4Stage6thOrderLobattoB_RKBT :
2058  virtual public RKButcherTableauDefaultBase<Scalar>
2059 {
2060  public:
2061  Implicit4Stage6thOrderLobattoB_RKBT()
2062  {
2063 
2064  std::ostringstream myDescription;
2065  myDescription << Implicit4Stage6thOrderLobattoB_name() << "\n"
2066  << "A-stable\n"
2067  << "Solving Ordinary Differential Equations II:\n"
2068  << "Stiff and Differential-Algebraic Problems,\n"
2069  << "2nd Revised Edition\n"
2070  << "E. Hairer and G. Wanner\n"
2071  << "Table 5.10, pg 76\n"
2072  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2073  << "A = [ 1/12 (-1-sqrt(5))/24 (-1+sqrt(5))/24 0 ]\n"
2074  << " [ 1/12 (25+sqrt(5))/120 (25-13*sqrt(5))/120 0 ]\n"
2075  << " [ 1/12 (25+13*sqrt(5))/120 (25-sqrt(5))/120 0 ]\n"
2076  << " [ 1/12 (11-sqrt(5))/24 (11+sqrt(5))/24 0 ]\n"
2077  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2078  typedef ScalarTraits<Scalar> ST;
2079  int myNumStages = 4;
2080  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2081  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2082  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2083  Scalar zero = ST::zero();
2084  Scalar one = ST::one();
2085  myA(0,0) = as<Scalar>( one/(12*one) );
2086  myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
2087  myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
2088  myA(0,3) = zero;
2089  myA(1,0) = as<Scalar>( one/(12*one) );
2090  myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
2091  myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
2092  myA(1,3) = zero;
2093  myA(2,0) = as<Scalar>( one/(12*one) );
2094  myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
2095  myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
2096  myA(2,3) = zero;
2097  myA(3,0) = as<Scalar>( one/(12*one) );
2098  myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
2099  myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
2100  myA(3,3) = zero;
2101  myb(0) = as<Scalar>( one/(12*one) );
2102  myb(1) = as<Scalar>( 5*one/(12*one) );
2103  myb(2) = as<Scalar>( 5*one/(12*one) );
2104  myb(3) = as<Scalar>( one/(12*one) );
2105  myc(0) = zero;
2106  myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2107  myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2108  myc(3) = one;
2109  this->setMyDescription(myDescription.str());
2110  this->setMy_A(myA);
2111  this->setMy_b(myb);
2112  this->setMy_c(myc);
2113  this->setMy_order(6);
2114  }
2115 };
2116 
2117 
2118 template<class Scalar>
2119 class Implicit2Stage2ndOrderLobattoC_RKBT :
2120  virtual public RKButcherTableauDefaultBase<Scalar>
2121 {
2122  public:
2123  Implicit2Stage2ndOrderLobattoC_RKBT()
2124  {
2125 
2126  std::ostringstream myDescription;
2127  myDescription << Implicit2Stage2ndOrderLobattoC_name() << "\n"
2128  << "A-stable\n"
2129  << "Solving Ordinary Differential Equations II:\n"
2130  << "Stiff and Differential-Algebraic Problems,\n"
2131  << "2nd Revised Edition\n"
2132  << "E. Hairer and G. Wanner\n"
2133  << "Table 5.11, pg 76\n"
2134  << "c = [ 0 1 ]'\n"
2135  << "A = [ 1/2 -1/2 ]\n"
2136  << " [ 1/2 1/2 ]\n"
2137  << "b = [ 1/2 1/2 ]'" << std::endl;
2138  typedef ScalarTraits<Scalar> ST;
2139  int myNumStages = 2;
2140  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2141  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2142  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2143  Scalar zero = ST::zero();
2144  Scalar one = ST::one();
2145  myA(0,0) = as<Scalar>( one/(2*one) );
2146  myA(0,1) = as<Scalar>( -one/(2*one) );
2147  myA(1,0) = as<Scalar>( one/(2*one) );
2148  myA(1,1) = as<Scalar>( one/(2*one) );
2149  myb(0) = as<Scalar>( one/(2*one) );
2150  myb(1) = as<Scalar>( one/(2*one) );
2151  myc(0) = zero;
2152  myc(1) = one;
2153  this->setMyDescription(myDescription.str());
2154  this->setMy_A(myA);
2155  this->setMy_b(myb);
2156  this->setMy_c(myc);
2157  this->setMy_order(2);
2158  }
2159 };
2160 
2161 
2162 template<class Scalar>
2163 class Implicit3Stage4thOrderLobattoC_RKBT :
2164  virtual public RKButcherTableauDefaultBase<Scalar>
2165 {
2166  public:
2167  Implicit3Stage4thOrderLobattoC_RKBT()
2168  {
2169 
2170  std::ostringstream myDescription;
2171  myDescription << Implicit3Stage4thOrderLobattoC_name() << "\n"
2172  << "A-stable\n"
2173  << "Solving Ordinary Differential Equations II:\n"
2174  << "Stiff and Differential-Algebraic Problems,\n"
2175  << "2nd Revised Edition\n"
2176  << "E. Hairer and G. Wanner\n"
2177  << "Table 5.11, pg 76\n"
2178  << "c = [ 0 1/2 1 ]'\n"
2179  << "A = [ 1/6 -1/3 1/6 ]\n"
2180  << " [ 1/6 5/12 -1/12 ]\n"
2181  << " [ 1/6 2/3 1/6 ]\n"
2182  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2183  typedef ScalarTraits<Scalar> ST;
2184  int myNumStages = 3;
2185  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2186  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2187  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2188  Scalar zero = ST::zero();
2189  Scalar one = ST::one();
2190  myA(0,0) = as<Scalar>( one/(6*one) );
2191  myA(0,1) = as<Scalar>( -one/(3*one) );
2192  myA(0,2) = as<Scalar>( one/(6*one) );
2193  myA(1,0) = as<Scalar>( one/(6*one) );
2194  myA(1,1) = as<Scalar>( 5*one/(12*one) );
2195  myA(1,2) = as<Scalar>( -one/(12*one) );
2196  myA(2,0) = as<Scalar>( one/(6*one) );
2197  myA(2,1) = as<Scalar>( 2*one/(3*one) );
2198  myA(2,2) = as<Scalar>( one/(6*one) );
2199  myb(0) = as<Scalar>( one/(6*one) );
2200  myb(1) = as<Scalar>( 2*one/(3*one) );
2201  myb(2) = as<Scalar>( one/(6*one) );
2202  myc(0) = zero;
2203  myc(1) = as<Scalar>( one/(2*one) );
2204  myc(2) = one;
2205  this->setMyDescription(myDescription.str());
2206  this->setMy_A(myA);
2207  this->setMy_b(myb);
2208  this->setMy_c(myc);
2209  this->setMy_order(4);
2210  }
2211 };
2212 
2213 
2214 template<class Scalar>
2215 class Implicit4Stage6thOrderLobattoC_RKBT :
2216  virtual public RKButcherTableauDefaultBase<Scalar>
2217 {
2218  public:
2219  Implicit4Stage6thOrderLobattoC_RKBT()
2220  {
2221 
2222  std::ostringstream myDescription;
2223  myDescription << Implicit4Stage6thOrderLobattoC_name() << "\n"
2224  << "A-stable\n"
2225  << "Solving Ordinary Differential Equations II:\n"
2226  << "Stiff and Differential-Algebraic Problems,\n"
2227  << "2nd Revised Edition\n"
2228  << "E. Hairer and G. Wanner\n"
2229  << "Table 5.12, pg 76\n"
2230  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2231  << "A = [ 1/12 -sqrt(5)/12 sqrt(5)/12 -1/12 ]\n"
2232  << " [ 1/12 1/4 (10-7*sqrt(5))/60 sqrt(5)/60 ]\n"
2233  << " [ 1/12 (10+7*sqrt(5))/60 1/4 -sqrt(5)/60 ]\n"
2234  << " [ 1/12 5/12 5/12 1/12 ]\n"
2235  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2236  typedef ScalarTraits<Scalar> ST;
2237  int myNumStages = 4;
2238  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2239  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2240  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2241  Scalar zero = ST::zero();
2242  Scalar one = ST::one();
2243  myA(0,0) = as<Scalar>( one/(12*one) );
2244  myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
2245  myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
2246  myA(0,3) = as<Scalar>( -one/(12*one) );
2247  myA(1,0) = as<Scalar>( one/(12*one) );
2248  myA(1,1) = as<Scalar>( one/(4*one) );
2249  myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
2250  myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
2251  myA(2,0) = as<Scalar>( one/(12*one) );
2252  myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
2253  myA(2,2) = as<Scalar>( one/(4*one) );
2254  myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
2255  myA(3,0) = as<Scalar>( one/(12*one) );
2256  myA(3,1) = as<Scalar>( 5*one/(12*one) );
2257  myA(3,2) = as<Scalar>( 5*one/(12*one) );
2258  myA(3,3) = as<Scalar>( one/(12*one) );
2259  myb(0) = as<Scalar>( one/(12*one) );
2260  myb(1) = as<Scalar>( 5*one/(12*one) );
2261  myb(2) = as<Scalar>( 5*one/(12*one) );
2262  myb(3) = as<Scalar>( one/(12*one) );
2263  myc(0) = zero;
2264  myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2265  myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2266  myc(3) = one;
2267  this->setMyDescription(myDescription.str());
2268  this->setMy_A(myA);
2269  this->setMy_b(myb);
2270  this->setMy_c(myc);
2271  this->setMy_order(6);
2272  }
2273 };
2274 
2275 
2276 
2277 template<class Scalar>
2278 class SDIRK5Stage5thOrder_RKBT :
2279  virtual public RKButcherTableauDefaultBase<Scalar>
2280 {
2281  public:
2282  SDIRK5Stage5thOrder_RKBT()
2283  {
2284 
2285  std::ostringstream myDescription;
2286  myDescription << SDIRK5Stage5thOrder_name() << "\n"
2287  << "A-stable\n"
2288  << "Solving Ordinary Differential Equations II:\n"
2289  << "Stiff and Differential-Algebraic Problems,\n"
2290  << "2nd Revised Edition\n"
2291  << "E. Hairer and G. Wanner\n"
2292  << "pg101 \n"
2293  << "c = [ (6-sqrt(6))/10 ]\n"
2294  << " [ (6+9*sqrt(6))/35 ]\n"
2295  << " [ 1 ]\n"
2296  << " [ (4-sqrt(6))/10 ]\n"
2297  << " [ (4+sqrt(6))/10 ]\n"
2298  << "A = [ A1 A2 A3 A4 A5 ]\n"
2299  << " A1 = [ (6-sqrt(6))/10 ]\n"
2300  << " [ (-6+5*sqrt(6))/14 ]\n"
2301  << " [ (888+607*sqrt(6))/2850 ]\n"
2302  << " [ (3153-3082*sqrt(6))/14250 ]\n"
2303  << " [ (-32583+14638*sqrt(6))/71250 ]\n"
2304  << " A2 = [ 0 ]\n"
2305  << " [ (6-sqrt(6))/10 ]\n"
2306  << " [ (126-161*sqrt(6))/1425 ]\n"
2307  << " [ (3213+1148*sqrt(6))/28500 ]\n"
2308  << " [ (-17199+364*sqrt(6))/142500 ]\n"
2309  << " A3 = [ 0 ]\n"
2310  << " [ 0 ]\n"
2311  << " [ (6-sqrt(6))/10 ]\n"
2312  << " [ (-267+88*sqrt(6))/500 ]\n"
2313  << " [ (1329-544*sqrt(6))/2500 ]\n"
2314  << " A4 = [ 0 ]\n"
2315  << " [ 0 ]\n"
2316  << " [ 0 ]\n"
2317  << " [ (6-sqrt(6))/10 ]\n"
2318  << " [ (-96+131*sqrt(6))/625 ]\n"
2319  << " A5 = [ 0 ]\n"
2320  << " [ 0 ]\n"
2321  << " [ 0 ]\n"
2322  << " [ 0 ]\n"
2323  << " [ (6-sqrt(6))/10 ]\n"
2324  << "b = [ 0 ]\n"
2325  << " [ 0 ]\n"
2326  << " [ 1/9 ]\n"
2327  << " [ (16-sqrt(6))/36 ]\n"
2328  << " [ (16+sqrt(6))/36 ]" << std::endl;
2329  typedef ScalarTraits<Scalar> ST;
2330  int myNumStages = 5;
2331  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2332  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2333  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2334  Scalar zero = ST::zero();
2335  Scalar one = ST::one();
2336  Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
2337  Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
2338  myA(0,0) = gamma;
2339  myA(0,1) = zero;
2340  myA(0,2) = zero;
2341  myA(0,3) = zero;
2342  myA(0,4) = zero;
2343 
2344  myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
2345  myA(1,1) = gamma;
2346  myA(1,2) = zero;
2347  myA(1,3) = zero;
2348  myA(1,4) = zero;
2349 
2350  myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
2351  myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
2352  myA(2,2) = gamma;
2353  myA(2,3) = zero;
2354  myA(2,4) = zero;
2355 
2356  myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
2357  myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
2358  myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
2359  myA(3,3) = gamma;
2360  myA(3,4) = zero;
2361 
2362  myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
2363  myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
2364  myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
2365  myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
2366  myA(4,4) = gamma;
2367 
2368  myb(0) = zero;
2369  myb(1) = zero;
2370  myb(2) = as<Scalar>( one/(9*one) );
2371  myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
2372  myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
2373 
2374  myc(0) = gamma;
2375  myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
2376  myc(2) = one;
2377  myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
2378  myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
2379 
2380  this->setMyDescription(myDescription.str());
2381  this->setMy_A(myA);
2382  this->setMy_b(myb);
2383  this->setMy_c(myc);
2384  this->setMy_order(5);
2385  }
2386 };
2387 
2388 
2389 template<class Scalar>
2390 class SDIRK5Stage4thOrder_RKBT :
2391  virtual public RKButcherTableauDefaultBase<Scalar>
2392 {
2393  public:
2394  SDIRK5Stage4thOrder_RKBT()
2395  {
2396 
2397  std::ostringstream myDescription;
2398  myDescription << SDIRK5Stage4thOrder_name() << "\n"
2399  << "L-stable\n"
2400  << "Solving Ordinary Differential Equations II:\n"
2401  << "Stiff and Differential-Algebraic Problems,\n"
2402  << "2nd Revised Edition\n"
2403  << "E. Hairer and G. Wanner\n"
2404  << "pg100 \n"
2405  << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
2406  << "A = [ 1/4 ]\n"
2407  << " [ 1/2 1/4 ]\n"
2408  << " [ 17/50 -1/25 1/4 ]\n"
2409  << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
2410  << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
2411  << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'\n"
2412  << "b' = [ 59/48 -17/96 225/32 -85/12 0 ]'" << std::endl;
2413  typedef ScalarTraits<Scalar> ST;
2414  int myNumStages = 5;
2415  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2416  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2417  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2418  Scalar zero = ST::zero();
2419  Scalar one = ST::one();
2420  Scalar onequarter = as<Scalar>( one/(4*one) );
2421  myA(0,0) = onequarter;
2422  myA(0,1) = zero;
2423  myA(0,2) = zero;
2424  myA(0,3) = zero;
2425  myA(0,4) = zero;
2426 
2427  myA(1,0) = as<Scalar>( one / (2*one) );
2428  myA(1,1) = onequarter;
2429  myA(1,2) = zero;
2430  myA(1,3) = zero;
2431  myA(1,4) = zero;
2432 
2433  myA(2,0) = as<Scalar>( 17*one/(50*one) );
2434  myA(2,1) = as<Scalar>( -one/(25*one) );
2435  myA(2,2) = onequarter;
2436  myA(2,3) = zero;
2437  myA(2,4) = zero;
2438 
2439  myA(3,0) = as<Scalar>( 371*one/(1360*one) );
2440  myA(3,1) = as<Scalar>( -137*one/(2720*one) );
2441  myA(3,2) = as<Scalar>( 15*one/(544*one) );
2442  myA(3,3) = onequarter;
2443  myA(3,4) = zero;
2444 
2445  myA(4,0) = as<Scalar>( 25*one/(24*one) );
2446  myA(4,1) = as<Scalar>( -49*one/(48*one) );
2447  myA(4,2) = as<Scalar>( 125*one/(16*one) );
2448  myA(4,3) = as<Scalar>( -85*one/(12*one) );
2449  myA(4,4) = onequarter;
2450 
2451  myb(0) = as<Scalar>( 25*one/(24*one) );
2452  myb(1) = as<Scalar>( -49*one/(48*one) );
2453  myb(2) = as<Scalar>( 125*one/(16*one) );
2454  myb(3) = as<Scalar>( -85*one/(12*one) );
2455  myb(4) = onequarter;
2456 
2457  /*
2458  // Alternate version
2459  myb(0) = as<Scalar>( 59*one/(48*one) );
2460  myb(1) = as<Scalar>( -17*one/(96*one) );
2461  myb(2) = as<Scalar>( 225*one/(32*one) );
2462  myb(3) = as<Scalar>( -85*one/(12*one) );
2463  myb(4) = zero;
2464  */
2465  myc(0) = onequarter;
2466  myc(1) = as<Scalar>( 3*one/(4*one) );
2467  myc(2) = as<Scalar>( 11*one/(20*one) );
2468  myc(3) = as<Scalar>( one/(2*one) );
2469  myc(4) = one;
2470 
2471  this->setMyDescription(myDescription.str());
2472  this->setMy_A(myA);
2473  this->setMy_b(myb);
2474  this->setMy_c(myc);
2475  this->setMy_order(4);
2476  }
2477 };
2478 
2479 
2480 template<class Scalar>
2481 class SDIRK3Stage4thOrder_RKBT :
2482  virtual public RKButcherTableauDefaultBase<Scalar>
2483 {
2484  public:
2485  SDIRK3Stage4thOrder_RKBT()
2486  {
2487 
2488  std::ostringstream myDescription;
2489  myDescription << SDIRK3Stage4thOrder_name() << "\n"
2490  << "A-stable\n"
2491  << "Solving Ordinary Differential Equations II:\n"
2492  << "Stiff and Differential-Algebraic Problems,\n"
2493  << "2nd Revised Edition\n"
2494  << "E. Hairer and G. Wanner\n"
2495  << "pg100 \n"
2496  << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
2497  << "delta = 1/(6*(2*gamma-1)^2)\n"
2498  << "c = [ gamma 1/2 1-gamma ]'\n"
2499  << "A = [ gamma ]\n"
2500  << " [ 1/2-gamma gamma ]\n"
2501  << " [ 2*gamma 1-4*gamma gamma ]\n"
2502  << "b = [ delta 1-2*delta delta ]'" << std::endl;
2503  typedef ScalarTraits<Scalar> ST;
2504  int myNumStages = 3;
2505  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2506  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2507  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2508  Scalar zero = ST::zero();
2509  Scalar one = ST::one();
2510  Scalar pi = as<Scalar>(4*one)*std::atan(one);
2511  Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
2512  Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
2513  myA(0,0) = gamma;
2514  myA(0,1) = zero;
2515  myA(0,2) = zero;
2516 
2517  myA(1,0) = as<Scalar>( one/(2*one) - gamma );
2518  myA(1,1) = gamma;
2519  myA(1,2) = zero;
2520 
2521  myA(2,0) = as<Scalar>( 2*gamma );
2522  myA(2,1) = as<Scalar>( one - 4*gamma );
2523  myA(2,2) = gamma;
2524 
2525  myb(0) = delta;
2526  myb(1) = as<Scalar>( one-2*delta );
2527  myb(2) = delta;
2528 
2529  myc(0) = gamma;
2530  myc(1) = as<Scalar>( one/(2*one) );
2531  myc(2) = as<Scalar>( one - gamma );
2532 
2533  this->setMyDescription(myDescription.str());
2534  this->setMy_A(myA);
2535  this->setMy_b(myb);
2536  this->setMy_c(myc);
2537  this->setMy_order(4);
2538  }
2539 };
2540 
2541 
2542 } // namespace Rythmos
2543 
2544 
2545 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP