30 #ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HPP 31 #define RYTHMOS_RK_BUTCHER_TABLEAU_HPP 35 #pragma clang system_header 38 #include "Rythmos_Types.hpp" 39 #include "Rythmos_RKButcherTableauBase.hpp" 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" 50 #include "Thyra_ProductVectorBase.hpp" 56 inline const std::string RKBT_ForwardEuler_name() {
return "Forward Euler"; }
57 inline const std::string RKBT_BackwardEuler_name() {
return "Backward Euler"; }
58 inline const std::string Explicit4Stage_name() {
return "Explicit 4 Stage"; }
59 inline const std::string Explicit3_8Rule_name() {
return "Explicit 3/8 Rule"; }
61 inline const std::string ExplicitTrapezoidal_name() {
return "Explicit Trapezoidal"; }
62 inline const std::string Explicit2Stage2ndOrderRunge_name() {
return "Explicit 2 Stage 2nd order by Runge"; }
63 inline const std::string Explicit2Stage2ndOrderTVD_name() {
return "Explicit 2 Stage 2nd order TVD"; }
64 inline const std::string Explicit3Stage3rdOrderHeun_name() {
return "Explicit 3 Stage 3rd order by Heun"; }
65 inline const std::string Explicit3Stage3rdOrder_name() {
return "Explicit 3 Stage 3rd order"; }
66 inline const std::string Explicit3Stage3rdOrderTVD_name() {
return "Explicit 3 Stage 3rd order TVD"; }
67 inline const std::string Explicit4Stage3rdOrderRunge_name() {
return "Explicit 4 Stage 3rd order by Runge"; }
68 inline const std::string Explicit5Stage3rdOrderKandG_name() {
return "Explicit 5 Stage 3rd order by Kinnmark and Gray"; }
70 inline const std::string IRK1StageTheta_name() {
return "IRK 1 Stage Theta Method"; }
71 inline const std::string IRK2StageTheta_name() {
return "IRK 2 Stage Theta Method"; }
72 inline const std::string Implicit1Stage2ndOrderGauss_name() {
return "Implicit 1 Stage 2nd order Gauss"; }
73 inline const std::string Implicit2Stage4thOrderGauss_name() {
return "Implicit 2 Stage 4th order Gauss"; }
74 inline const std::string Implicit3Stage6thOrderGauss_name() {
return "Implicit 3 Stage 6th order Gauss"; }
76 inline const std::string Implicit1Stage1stOrderRadauA_name() {
return "Implicit 1 Stage 1st order Radau left"; }
77 inline const std::string Implicit2Stage3rdOrderRadauA_name() {
return "Implicit 2 Stage 3rd order Radau left"; }
78 inline const std::string Implicit3Stage5thOrderRadauA_name() {
return "Implicit 3 Stage 5th order Radau left"; }
80 inline const std::string Implicit1Stage1stOrderRadauB_name() {
return "Implicit 1 Stage 1st order Radau right"; }
81 inline const std::string Implicit2Stage3rdOrderRadauB_name() {
return "Implicit 2 Stage 3rd order Radau right"; }
82 inline const std::string Implicit3Stage5thOrderRadauB_name() {
return "Implicit 3 Stage 5th order Radau right"; }
84 inline const std::string Implicit2Stage2ndOrderLobattoA_name() {
return "Implicit 2 Stage 2nd order Lobatto A"; }
85 inline const std::string Implicit3Stage4thOrderLobattoA_name() {
return "Implicit 3 Stage 4th order Lobatto A"; }
86 inline const std::string Implicit4Stage6thOrderLobattoA_name() {
return "Implicit 4 Stage 6th order Lobatto A"; }
88 inline const std::string Implicit2Stage2ndOrderLobattoB_name() {
return "Implicit 2 Stage 2nd order Lobatto B"; }
89 inline const std::string Implicit3Stage4thOrderLobattoB_name() {
return "Implicit 3 Stage 4th order Lobatto B"; }
90 inline const std::string Implicit4Stage6thOrderLobattoB_name() {
return "Implicit 4 Stage 6th order Lobatto B"; }
92 inline const std::string Implicit2Stage2ndOrderLobattoC_name() {
return "Implicit 2 Stage 2nd order Lobatto C"; }
93 inline const std::string Implicit3Stage4thOrderLobattoC_name() {
return "Implicit 3 Stage 4th order Lobatto C"; }
94 inline const std::string Implicit4Stage6thOrderLobattoC_name() {
return "Implicit 4 Stage 6th order Lobatto C"; }
96 inline const std::string Implicit2Stage4thOrderHammerHollingsworth_name() {
return "Implicit 2 Stage 4th Order Hammer & Hollingsworth"; }
97 inline const std::string Implicit3Stage6thOrderKuntzmannButcher_name() {
return "Implicit 3 Stage 6th Order Kuntzmann & Butcher"; }
98 inline const std::string Implicit4Stage8thOrderKuntzmannButcher_name() {
return "Implicit 4 Stage 8th Order Kuntzmann & Butcher"; }
100 inline const std::string DIRK2Stage3rdOrder_name() {
return "Diagonal IRK 2 Stage 3rd order"; }
102 inline const std::string SDIRK2Stage2ndOrder_name() {
return "Singly Diagonal IRK 2 Stage 2nd order"; }
103 inline const std::string SDIRK2Stage3rdOrder_name() {
return "Singly Diagonal IRK 2 Stage 3rd order"; }
104 inline const std::string SDIRK5Stage5thOrder_name() {
return "Singly Diagonal IRK 5 Stage 5th order"; }
105 inline const std::string SDIRK5Stage4thOrder_name() {
return "Singly Diagonal IRK 5 Stage 4th order"; }
106 inline const std::string SDIRK3Stage4thOrder_name() {
return "Singly Diagonal IRK 3 Stage 4th order"; }
108 template<
class Scalar>
109 class RKButcherTableauDefaultBase :
110 virtual public RKButcherTableauBase<Scalar>,
111 virtual public Teuchos::ParameterListAcceptorDefaultBase
115 virtual int numStages()
const {
return A_.numRows(); }
117 virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A()
const {
return A_; }
119 virtual const Teuchos::SerialDenseVector<int,Scalar>& b()
const {
return b_; }
121 virtual const Teuchos::SerialDenseVector<int,Scalar>& bhat()
const {
return bhat_ ; }
123 virtual const Teuchos::SerialDenseVector<int,Scalar>& c()
const {
return c_; }
125 virtual int order()
const {
return order_; }
127 virtual bool isEmbeddedMethod()
const {
return isEmbedded_; }
129 virtual void setDescription(std::string longDescription) { longDescription_ = longDescription; }
132 virtual void initialize(
133 const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
134 const Teuchos::SerialDenseVector<int,Scalar>& b_in,
135 const Teuchos::SerialDenseVector<int,Scalar>& c_in,
137 const std::string& longDescription_in,
138 bool isEmbedded =
false,
139 const Teuchos::SerialDenseVector<int,Scalar>& bhat_in = Teuchos::SerialDenseVector<int,Scalar>()
142 const int numStages_in = A_in.numRows();
143 TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
144 TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
145 TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
146 TEUCHOS_ASSERT_EQUALITY( c_in.length(), numStages_in );
147 TEUCHOS_ASSERT( order_in > 0 );
152 longDescription_ = longDescription_in;
156 TEUCHOS_ASSERT_EQUALITY( bhat_in.length(), numStages_in );
166 virtual void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
168 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
169 paramList->validateParameters(*this->getValidParameters());
170 Teuchos::readVerboseObjectSublist(&*paramList,
this);
171 setMyParamList(paramList);
175 virtual RCP<const Teuchos::ParameterList> getValidParameters()
const 177 if (is_null(validPL_)) {
178 validPL_ = Teuchos::parameterList();
179 validPL_->set(
"Description",
"",this->getMyDescription());
180 Teuchos::setupVerboseObjectSublist(&*validPL_);
191 virtual std::string description()
const {
return "Rythmos::RKButcherTableauDefaultBase"; }
194 virtual void describe(
195 Teuchos::FancyOStream &out,
196 const Teuchos::EVerbosityLevel verbLevel
199 if (verbLevel != Teuchos::VERB_NONE) {
200 out << this->description() << std::endl;
201 out << this->getMyDescription() << std::endl;
202 out <<
"number of Stages = " << this->numStages() << std::endl;
203 out <<
"A = " << this->A() << std::endl;
204 out <<
"b = " << this->b() << std::endl;
205 out <<
"c = " << this->c() << std::endl;
206 out <<
"order = " << this->order() << std::endl;
213 void setMyDescription(std::string longDescription) { longDescription_ = longDescription; }
214 const std::string& getMyDescription()
const {
return longDescription_; }
216 void setMy_A(
const Teuchos::SerialDenseMatrix<int,Scalar>& new_A) { A_ = new_A; }
217 void setMy_b(
const Teuchos::SerialDenseVector<int,Scalar>& new_b) { b_ = new_b; }
218 void setMy_c(
const Teuchos::SerialDenseVector<int,Scalar>& new_c) { c_ = new_c; }
219 void setMy_order(
const int& new_order) { order_ = new_order; }
221 void setMyValidParameterList(
const RCP<ParameterList> validPL ) { validPL_ = validPL; }
222 RCP<ParameterList> getMyNonconstValidParameterList() {
return validPL_; }
225 Teuchos::SerialDenseMatrix<int,Scalar> A_;
226 Teuchos::SerialDenseVector<int,Scalar> b_;
227 Teuchos::SerialDenseVector<int,Scalar> c_;
229 std::string longDescription_;
230 mutable RCP<ParameterList> validPL_;
233 Teuchos::SerialDenseVector<int,Scalar> bhat_;
234 bool isEmbedded_ =
false;
239 template<
class Scalar>
240 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
242 return(rcp(
new RKButcherTableauDefaultBase<Scalar>()));
246 template<
class Scalar>
247 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau(
248 const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
249 const Teuchos::SerialDenseVector<int,Scalar>& b_in,
250 const Teuchos::SerialDenseVector<int,Scalar>& c_in,
252 const std::string& description_in =
"" 255 RCP<RKButcherTableauDefaultBase<Scalar> > rkbt = rcp(
new RKButcherTableauDefaultBase<Scalar>());
256 rkbt->initialize(A_in,b_in,c_in,order_in,description_in);
261 template<
class Scalar>
262 class BackwardEuler_RKBT :
263 virtual public RKButcherTableauDefaultBase<Scalar>
268 std::ostringstream myDescription;
269 myDescription << RKBT_BackwardEuler_name() <<
"\n" 272 <<
"b = [ 1 ]'" << std::endl;
273 typedef ScalarTraits<Scalar> ST;
274 Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
275 myA(0,0) = ST::one();
276 Teuchos::SerialDenseVector<int,Scalar> myb(1);
278 Teuchos::SerialDenseVector<int,Scalar> myc(1);
281 this->setMyDescription(myDescription.str());
285 this->setMy_order(1);
291 template<
class Scalar>
292 class ForwardEuler_RKBT :
293 virtual public RKButcherTableauDefaultBase<Scalar>
299 std::ostringstream myDescription;
300 myDescription << RKBT_ForwardEuler_name() <<
"\n" 303 <<
"b = [ 1 ]'" << std::endl;
304 typedef ScalarTraits<Scalar> ST;
305 Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
306 Teuchos::SerialDenseVector<int,Scalar> myb(1);
308 Teuchos::SerialDenseVector<int,Scalar> myc(1);
310 this->setMyDescription(myDescription.str());
314 this->setMy_order(1);
319 template<
class Scalar>
320 class Explicit4Stage4thOrder_RKBT :
321 virtual public RKButcherTableauDefaultBase<Scalar>
324 Explicit4Stage4thOrder_RKBT()
326 std::ostringstream myDescription;
327 myDescription << Explicit4Stage_name() <<
"\n" 328 <<
"\"The\" Runge-Kutta Method (explicit):\n" 329 <<
"Solving Ordinary Differential Equations I:\n" 330 <<
"Nonstiff Problems, 2nd Revised Edition\n" 331 <<
"E. Hairer, S.P. Norsett, G. Wanner\n" 332 <<
"Table 1.2, pg 138\n" 333 <<
"c = [ 0 1/2 1/2 1 ]'\n" 338 <<
"b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
339 typedef ScalarTraits<Scalar> ST;
340 Scalar one = ST::one();
341 Scalar zero = ST::zero();
342 Scalar onehalf = ST::one()/(2*ST::one());
343 Scalar onesixth = ST::one()/(6*ST::one());
344 Scalar onethird = ST::one()/(3*ST::one());
347 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
348 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
349 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
384 this->setMyDescription(myDescription.str());
388 this->setMy_order(4);
393 template<
class Scalar>
394 class Explicit3_8Rule_RKBT :
395 virtual public RKButcherTableauDefaultBase<Scalar>
398 Explicit3_8Rule_RKBT()
401 std::ostringstream myDescription;
402 myDescription << Explicit3_8Rule_name() <<
"\n" 403 <<
"Solving Ordinary Differential Equations I:\n" 404 <<
"Nonstiff Problems, 2nd Revised Edition\n" 405 <<
"E. Hairer, S.P. Norsett, G. Wanner\n" 406 <<
"Table 1.2, pg 138\n" 407 <<
"c = [ 0 1/3 2/3 1 ]'\n" 412 <<
"b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
413 typedef ScalarTraits<Scalar> ST;
415 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
416 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
417 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
419 Scalar one = ST::one();
420 Scalar zero = ST::zero();
421 Scalar one_third = as<Scalar>(ST::one()/(3*ST::one()));
422 Scalar two_third = as<Scalar>(2*ST::one()/(3*ST::one()));
423 Scalar one_eighth = as<Scalar>(ST::one()/(8*ST::one()));
424 Scalar three_eighth = as<Scalar>(3*ST::one()/(8*ST::one()));
432 myA(1,0) = one_third;
437 myA(2,0) = as<Scalar>(-one_third);
443 myA(3,1) = as<Scalar>(-one);
449 myb(1) = three_eighth;
450 myb(2) = three_eighth;
459 this->setMyDescription(myDescription.str());
463 this->setMy_order(4);
468 template<
class Scalar>
469 class Explicit4Stage3rdOrderRunge_RKBT :
470 virtual public RKButcherTableauDefaultBase<Scalar>
473 Explicit4Stage3rdOrderRunge_RKBT()
476 std::ostringstream myDescription;
477 myDescription << Explicit4Stage3rdOrderRunge_name() <<
"\n" 478 <<
"Solving Ordinary Differential Equations I:\n" 479 <<
"Nonstiff Problems, 2nd Revised Edition\n" 480 <<
"E. Hairer, S.P. Norsett, G. Wanner\n" 481 <<
"Table 1.1, pg 135\n" 482 <<
"c = [ 0 1/2 1 1 ]'\n" 487 <<
"b = [ 1/6 2/3 0 1/6 ]'" << std::endl;
488 typedef ScalarTraits<Scalar> ST;
490 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
491 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
492 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
494 Scalar one = ST::one();
495 Scalar onehalf = ST::one()/(2*ST::one());
496 Scalar onesixth = ST::one()/(6*ST::one());
497 Scalar twothirds = 2*ST::one()/(3*ST::one());
498 Scalar zero = ST::zero();
533 this->setMyDescription(myDescription.str());
537 this->setMy_order(3);
541 template<
class Scalar>
542 class Explicit5Stage3rdOrderKandG_RKBT :
543 virtual public RKButcherTableauDefaultBase<Scalar>
546 Explicit5Stage3rdOrderKandG_RKBT()
549 std::ostringstream myDescription;
550 myDescription << Explicit5Stage3rdOrderKandG_name() <<
"\n" 551 <<
"Kinnmark & Gray 5 stage, 3rd order scheme \n" 552 <<
"Modified by P. Ullrich. From the prim_advance_mod.F90 \n" 553 <<
"routine in the HOMME atmosphere model code.\n" 554 <<
"c = [ 0 1/5 1/5 1/3 2/3 ]'\n" 558 <<
" [ 0 0 1/3 0 ]\n" 559 <<
" [ 0 0 0 2/3 0 ]\n" 560 <<
"b = [ 1/4 0 0 0 3/4 ]'" << std::endl;
561 typedef ScalarTraits<Scalar> ST;
563 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
564 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
565 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
567 Scalar one = ST::one();
568 Scalar onefifth = ST::one()/(5*ST::one());
569 Scalar onefourth = ST::one()/(4*ST::one());
570 Scalar onethird = ST::one()/(3*ST::one());
571 Scalar twothirds = 2*ST::one()/(3*ST::one());
572 Scalar threefourths = 3*ST::one()/(4*ST::one());
573 Scalar zero = ST::zero();
603 myA(4,3) = twothirds;
611 myb(4) = threefourths;
620 this->setMyDescription(myDescription.str());
624 this->setMy_order(3);
629 template<
class Scalar>
630 class Explicit3Stage3rdOrder_RKBT :
631 virtual public RKButcherTableauDefaultBase<Scalar>
634 Explicit3Stage3rdOrder_RKBT()
637 std::ostringstream myDescription;
638 myDescription << Explicit3Stage3rdOrder_name() <<
"\n" 639 <<
"c = [ 0 1/2 1 ]'\n" 643 <<
"b = [ 1/6 4/6 1/6 ]'" << std::endl;
644 typedef ScalarTraits<Scalar> ST;
645 Scalar one = ST::one();
646 Scalar two = as<Scalar>(2*ST::one());
647 Scalar zero = ST::zero();
648 Scalar onehalf = ST::one()/(2*ST::one());
649 Scalar onesixth = ST::one()/(6*ST::one());
650 Scalar foursixth = 4*ST::one()/(6*ST::one());
653 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
654 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
655 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
680 this->setMyDescription(myDescription.str());
684 this->setMy_order(3);
688 template<
class Scalar>
689 class Explicit2Stage2ndOrderTVD_RKBT :
690 virtual public RKButcherTableauDefaultBase<Scalar>
693 Explicit2Stage2ndOrderTVD_RKBT()
696 std::ostringstream myDescription;
697 myDescription << Explicit2Stage2ndOrderTVD_name() <<
"\n" 698 <<
"Sigal Gottlieb, David Ketcheson and Chi-Wang Shu\n" 699 <<
"`Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations'\n" 700 <<
"World Scientific, 2011\n" 705 <<
"b = [ 1/2 1/2]'\n" 706 <<
"This is also written in the following set of updates.\n" 707 <<
"u1 = u^n + dt L(u^n)\n" 708 <<
"u^(n+1) = u^n/2 + u1/2 + dt L(u1)/2" 710 typedef ScalarTraits<Scalar> ST;
711 Scalar one = ST::one();
712 Scalar zero = ST::zero();
713 Scalar onehalf = ST::one()/(2*ST::one());
716 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
717 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
718 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
735 this->setMyDescription(myDescription.str());
739 this->setMy_order(2);
743 template<
class Scalar>
744 class Explicit3Stage3rdOrderTVD_RKBT :
745 virtual public RKButcherTableauDefaultBase<Scalar>
748 Explicit3Stage3rdOrderTVD_RKBT()
751 std::ostringstream myDescription;
752 myDescription << Explicit3Stage3rdOrderTVD_name() <<
"\n" 753 <<
"Sigal Gottlieb and Chi-Wang Shu\n" 754 <<
"`Total Variation Diminishing Runge-Kutta Schemes'\n" 755 <<
"Mathematics of Computation\n" 756 <<
"Volume 67, Number 221, January 1998, pp. 73-85\n" 757 <<
"c = [ 0 1 1/2 ]'\n" 760 <<
" [ 1/4 1/4 0 ]\n" 761 <<
"b = [ 1/6 1/6 4/6 ]'\n" 762 <<
"This is also written in the following set of updates.\n" 763 <<
"u1 = u^n + dt L(u^n)\n" 764 <<
"u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n" 765 <<
"u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3" 767 typedef ScalarTraits<Scalar> ST;
768 Scalar one = ST::one();
769 Scalar zero = ST::zero();
770 Scalar onehalf = ST::one()/(2*ST::one());
771 Scalar onefourth = ST::one()/(4*ST::one());
772 Scalar onesixth = ST::one()/(6*ST::one());
773 Scalar foursixth = 4*ST::one()/(6*ST::one());
776 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
777 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
778 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
789 myA(2,0) = onefourth;
790 myA(2,1) = onefourth;
803 this->setMyDescription(myDescription.str());
807 this->setMy_order(3);
812 template<
class Scalar>
813 class Explicit3Stage3rdOrderHeun_RKBT :
814 virtual public RKButcherTableauDefaultBase<Scalar>
817 Explicit3Stage3rdOrderHeun_RKBT()
819 std::ostringstream myDescription;
820 myDescription << Explicit3Stage3rdOrderHeun_name() <<
"\n" 821 <<
"Solving Ordinary Differential Equations I:\n" 822 <<
"Nonstiff Problems, 2nd Revised Edition\n" 823 <<
"E. Hairer, S.P. Norsett, G. Wanner\n" 824 <<
"Table 1.1, pg 135\n" 825 <<
"c = [ 0 1/3 2/3 ]'\n" 829 <<
"b = [ 1/4 0 3/4 ]'" << std::endl;
830 typedef ScalarTraits<Scalar> ST;
831 Scalar one = ST::one();
832 Scalar zero = ST::zero();
833 Scalar onethird = one/(3*one);
834 Scalar twothirds = 2*one/(3*one);
835 Scalar onefourth = one/(4*one);
836 Scalar threefourths = 3*one/(4*one);
839 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
840 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
841 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
853 myA(2,1) = twothirds;
859 myb(2) = threefourths;
866 this->setMyDescription(myDescription.str());
870 this->setMy_order(3);
875 template<
class Scalar>
876 class Explicit2Stage2ndOrderRunge_RKBT :
877 virtual public RKButcherTableauDefaultBase<Scalar>
880 Explicit2Stage2ndOrderRunge_RKBT()
882 std::ostringstream myDescription;
883 myDescription << Explicit2Stage2ndOrderRunge_name() <<
"\n" 884 <<
"Also known as Explicit Midpoint\n" 885 <<
"Solving Ordinary Differential Equations I:\n" 886 <<
"Nonstiff Problems, 2nd Revised Edition\n" 887 <<
"E. Hairer, S.P. Norsett, G. Wanner\n" 888 <<
"Table 1.1, pg 135\n" 889 <<
"c = [ 0 1/2 ]'\n" 892 <<
"b = [ 0 1 ]'" << std::endl;
893 typedef ScalarTraits<Scalar> ST;
894 Scalar one = ST::one();
895 Scalar zero = ST::zero();
896 Scalar onehalf = ST::one()/(2*ST::one());
899 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
900 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
901 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
918 this->setMyDescription(myDescription.str());
922 this->setMy_order(2);
927 template<
class Scalar>
928 class ExplicitTrapezoidal_RKBT :
929 virtual public RKButcherTableauDefaultBase<Scalar>
932 ExplicitTrapezoidal_RKBT()
934 std::ostringstream myDescription;
935 myDescription << ExplicitTrapezoidal_name() <<
"\n" 939 <<
"b = [ 1/2 1/2 ]'" << std::endl;
940 typedef ScalarTraits<Scalar> ST;
941 Scalar one = ST::one();
942 Scalar zero = ST::zero();
943 Scalar onehalf = ST::one()/(2*ST::one());
946 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
947 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
948 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
965 this->setMyDescription(myDescription.str());
969 this->setMy_order(2);
974 template<
class Scalar>
975 class SDIRK2Stage2ndOrder_RKBT :
976 virtual public RKButcherTableauDefaultBase<Scalar>
979 SDIRK2Stage2ndOrder_RKBT()
981 std::ostringstream myDescription;
982 myDescription << SDIRK2Stage2ndOrder_name() <<
"\n" 983 <<
"Computer Methods for ODEs and DAEs\n" 984 <<
"U. M. Ascher and L. R. Petzold\n" 986 <<
"gamma = (2+-sqrt(2))/2\n" 987 <<
"c = [ gamma 1 ]'\n" 988 <<
"A = [ gamma 0 ]\n" 989 <<
" [ 1-gamma gamma ]\n" 990 <<
"b = [ 1-gamma gamma ]'" << std::endl;
992 this->setMyDescription(myDescription.str());
993 typedef ScalarTraits<Scalar> ST;
994 Scalar one = ST::one();
995 gamma_default_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
996 gamma_ = gamma_default_;
999 RCP<ParameterList> validPL = Teuchos::parameterList();
1000 validPL->set(
"Description",
"",this->getMyDescription());
1001 validPL->set<
double>(
"gamma",gamma_default_,
1002 "The default value is gamma = (2-sqrt(2))/2. " 1003 "This will produce an L-stable 2nd order method with the stage " 1004 "times within the timestep. Other values of gamma will still " 1005 "produce an L-stable scheme, but will only be 1st order accurate.");
1006 Teuchos::setupVerboseObjectSublist(&*validPL);
1007 this->setMyValidParameterList(validPL);
1011 typedef ScalarTraits<Scalar> ST;
1012 int myNumStages = 2;
1013 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1014 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1015 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1016 Scalar one = ST::one();
1017 Scalar zero = ST::zero();
1020 myA(1,0) = as<Scalar>( one - gamma_ );
1022 myb(0) = as<Scalar>( one - gamma_ );
1030 this->setMy_order(2);
1032 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
1034 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1035 paramList->validateParameters(*this->getValidParameters());
1036 Teuchos::readVerboseObjectSublist(&*paramList,
this);
1037 gamma_ = paramList->get<
double>(
"gamma",gamma_default_);
1039 this->setMyParamList(paramList);
1042 Scalar gamma_default_;
1049 template<
class Scalar>
1050 class SDIRK2Stage3rdOrder_RKBT :
1051 virtual public RKButcherTableauDefaultBase<Scalar>
1054 SDIRK2Stage3rdOrder_RKBT()
1056 std::ostringstream myDescription;
1057 myDescription << SDIRK2Stage3rdOrder_name() <<
"\n" 1058 <<
"Solving Ordinary Differential Equations I:\n" 1059 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1060 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1061 <<
"Table 7.2, pg 207\n" 1062 <<
"gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n" 1063 <<
"gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n" 1064 <<
"c = [ gamma 1-gamma ]'\n" 1065 <<
"A = [ gamma 0 ]\n" 1066 <<
" [ 1-2*gamma gamma ]\n" 1067 <<
"b = [ 1/2 1/2 ]'" << std::endl;
1069 this->setMyDescription(myDescription.str());
1070 thirdOrderAStable_default_ =
true;
1071 secondOrderLStable_default_ =
false;
1072 thirdOrderAStable_ =
true;
1073 secondOrderLStable_ =
false;
1074 typedef ScalarTraits<Scalar> ST;
1075 Scalar one = ST::one();
1076 gamma_default_ = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1077 gamma_ = gamma_default_;
1080 RCP<ParameterList> validPL = Teuchos::parameterList();
1081 validPL->set(
"Description",
"",this->getMyDescription());
1082 validPL->set(
"3rd Order A-stable",thirdOrderAStable_default_,
1083 "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain " 1084 "a 3rd order A-stable scheme. '3rd Order A-stable' and " 1085 "'2nd Order L-stable' can not both be true.");
1086 validPL->set(
"2nd Order L-stable",secondOrderLStable_default_,
1087 "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain " 1088 "a 2nd order L-stable scheme. '3rd Order A-stable' and " 1089 "'2nd Order L-stable' can not both be true.");
1090 validPL->set(
"gamma",gamma_default_,
1091 "If both '3rd Order A-stable' and '2nd Order L-stable' " 1092 "are false, gamma will be used. The default value is the " 1093 "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
1095 Teuchos::setupVerboseObjectSublist(&*validPL);
1096 this->setMyValidParameterList(validPL);
1100 typedef ScalarTraits<Scalar> ST;
1101 int myNumStages = 2;
1102 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1103 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1104 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1105 Scalar one = ST::one();
1106 Scalar zero = ST::zero();
1108 if (thirdOrderAStable_)
1109 gamma = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1110 else if (secondOrderLStable_)
1111 gamma = as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
1116 myA(1,0) = as<Scalar>( one - 2*gamma );
1118 myb(0) = as<Scalar>( one/(2*one) );
1119 myb(1) = as<Scalar>( one/(2*one) );
1121 myc(1) = as<Scalar>( one - gamma );
1126 this->setMy_order(3);
1128 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
1130 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1131 paramList->validateParameters(*this->getValidParameters());
1132 Teuchos::readVerboseObjectSublist(&*paramList,
this);
1133 thirdOrderAStable_ = paramList->get(
"3rd Order A-stable",
1134 thirdOrderAStable_default_);
1135 secondOrderLStable_ = paramList->get(
"2nd Order L-stable",
1136 secondOrderLStable_default_);
1137 TEUCHOS_TEST_FOR_EXCEPTION(
1138 thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
1139 "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
1140 gamma_ = paramList->get(
"gamma",gamma_default_);
1143 this->setMyParamList(paramList);
1146 bool thirdOrderAStable_default_;
1147 bool thirdOrderAStable_;
1148 bool secondOrderLStable_default_;
1149 bool secondOrderLStable_;
1150 Scalar gamma_default_;
1155 template<
class Scalar>
1156 class DIRK2Stage3rdOrder_RKBT :
1157 virtual public RKButcherTableauDefaultBase<Scalar>
1160 DIRK2Stage3rdOrder_RKBT()
1163 std::ostringstream myDescription;
1164 myDescription << DIRK2Stage3rdOrder_name() <<
"\n" 1165 <<
"Hammer & Hollingsworth method\n" 1166 <<
"Solving Ordinary Differential Equations I:\n" 1167 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1168 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1169 <<
"Table 7.1, pg 205\n" 1170 <<
"c = [ 0 2/3 ]'\n" 1173 <<
"b = [ 1/4 3/4 ]'" << std::endl;
1174 typedef ScalarTraits<Scalar> ST;
1175 int myNumStages = 2;
1176 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1177 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1178 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1179 Scalar one = ST::one();
1180 Scalar zero = ST::zero();
1183 myA(1,0) = as<Scalar>( one/(3*one) );
1184 myA(1,1) = as<Scalar>( one/(3*one) );
1185 myb(0) = as<Scalar>( one/(4*one) );
1186 myb(1) = as<Scalar>( 3*one/(4*one) );
1188 myc(1) = as<Scalar>( 2*one/(3*one) );
1189 this->setMyDescription(myDescription.str());
1193 this->setMy_order(3);
1198 template<
class Scalar>
1199 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
1200 virtual public RKButcherTableauDefaultBase<Scalar>
1203 Implicit3Stage6thOrderKuntzmannButcher_RKBT()
1206 std::ostringstream myDescription;
1207 myDescription << Implicit3Stage6thOrderKuntzmannButcher_name() <<
"\n" 1208 <<
"Kuntzmann & Butcher method\n" 1209 <<
"Solving Ordinary Differential Equations I:\n" 1210 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1211 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1212 <<
"Table 7.4, pg 209\n" 1213 <<
"c = [ 1/2-sqrt(15)/10 1/2 1/2-sqrt(15)/10 ]'\n" 1214 <<
"A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n" 1215 <<
" [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n" 1216 <<
" [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n" 1217 <<
"b = [ 5/18 4/9 5/18 ]'" << std::endl;
1218 typedef ScalarTraits<Scalar> ST;
1219 int myNumStages = 3;
1220 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1221 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1222 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1223 Scalar one = ST::one();
1224 myA(0,0) = as<Scalar>( 5*one/(36*one) );
1225 myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
1226 myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
1227 myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
1228 myA(1,1) = as<Scalar>( 2*one/(9*one) );
1229 myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
1230 myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
1231 myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
1232 myA(2,2) = as<Scalar>( 5*one/(36*one) );
1233 myb(0) = as<Scalar>( 5*one/(18*one) );
1234 myb(1) = as<Scalar>( 4*one/(9*one) );
1235 myb(2) = as<Scalar>( 5*one/(18*one) );
1236 myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
1237 myc(1) = as<Scalar>( one/(2*one) );
1238 myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
1239 this->setMyDescription(myDescription.str());
1243 this->setMy_order(6);
1248 template<
class Scalar>
1249 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
1250 virtual public RKButcherTableauDefaultBase<Scalar>
1253 Implicit4Stage8thOrderKuntzmannButcher_RKBT()
1256 std::ostringstream myDescription;
1257 myDescription << Implicit4Stage8thOrderKuntzmannButcher_name() <<
"\n" 1258 <<
"Kuntzmann & Butcher method\n" 1259 <<
"Solving Ordinary Differential Equations I:\n" 1260 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1261 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1262 <<
"Table 7.5, pg 209\n" 1263 <<
"c = [ 1/2-w2 1/2-w2p 1/2+w2p 1/2+w2 ]'\n" 1264 <<
"A = [ w1 w1p-w3+w4p w1p-w3-w4p w1-w5 ]\n" 1265 <<
" [ w1-w3p+w4 w1p w1p-w5p w1-w3p-w4 ]\n" 1266 <<
" [ w1+w3p+w4 w1p+w5p w1p w1+w3p-w4 ]\n" 1267 <<
" [ w1+w5 w1p+w3+w4p w1p+w3-w4p w1 ]\n" 1268 <<
"b = [ 2*w1 2*w1p 2*w1p 2*w1 ]'\n" 1269 <<
"w1 = 1/8-sqrt(30)/144\n" 1270 <<
"w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n" 1271 <<
"w3 = w2*(1/6+sqrt(30)/24)\n" 1272 <<
"w4 = w2*(1/21+5*sqrt(30)/168)\n" 1274 <<
"w1p = 1/8+sqrt(30)/144\n" 1275 <<
"w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n" 1276 <<
"w3p = w2*(1/6-sqrt(30)/24)\n" 1277 <<
"w4p = w2*(1/21-5*sqrt(30)/168)\n" 1278 <<
"w5p = w2p-2*w3p" << std::endl;
1279 typedef ScalarTraits<Scalar> ST;
1280 int myNumStages = 4;
1281 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1282 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1283 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1284 Scalar one = ST::one();
1285 Scalar onehalf = as<Scalar>( one/(2*one) );
1286 Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
1287 Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
1288 Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
1289 Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
1290 Scalar w5 = as<Scalar>( w2-2*w3 );
1291 Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
1292 Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
1293 Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
1294 Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
1295 Scalar w5p = as<Scalar>( w2p-2*w3p );
1297 myA(0,1) = w1p-w3+w4p;
1298 myA(0,2) = w1p-w3-w4p;
1300 myA(1,0) = w1-w3p+w4;
1303 myA(1,3) = w1-w3p-w4;
1304 myA(2,0) = w1+w3p+w4;
1307 myA(2,3) = w1+w3p-w4;
1309 myA(3,1) = w1p+w3+w4p;
1310 myA(3,2) = w1p+w3-w4p;
1316 myc(0) = onehalf - w2;
1317 myc(1) = onehalf - w2p;
1318 myc(2) = onehalf + w2p;
1319 myc(3) = onehalf + w2;
1320 this->setMyDescription(myDescription.str());
1324 this->setMy_order(8);
1329 template<
class Scalar>
1330 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
1331 virtual public RKButcherTableauDefaultBase<Scalar>
1334 Implicit2Stage4thOrderHammerHollingsworth_RKBT()
1337 std::ostringstream myDescription;
1338 myDescription << Implicit2Stage4thOrderHammerHollingsworth_name() <<
"\n" 1339 <<
"Hammer & Hollingsworth method\n" 1340 <<
"Solving Ordinary Differential Equations I:\n" 1341 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1342 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1343 <<
"Table 7.3, pg 207\n" 1344 <<
"c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n" 1345 <<
"A = [ 1/4 1/4-sqrt(3)/6 ]\n" 1346 <<
" [ 1/4+sqrt(3)/6 1/4 ]\n" 1347 <<
"b = [ 1/2 1/2 ]'" << std::endl;
1348 typedef ScalarTraits<Scalar> ST;
1349 int myNumStages = 2;
1350 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1351 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1352 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1353 Scalar one = ST::one();
1354 Scalar onequarter = as<Scalar>( one/(4*one) );
1355 Scalar onehalf = as<Scalar>( one/(2*one) );
1356 myA(0,0) = onequarter;
1357 myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
1358 myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
1359 myA(1,1) = onequarter;
1362 myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
1363 myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
1364 this->setMyDescription(myDescription.str());
1368 this->setMy_order(4);
1373 template<
class Scalar>
1374 class IRK1StageTheta_RKBT :
1375 virtual public RKButcherTableauDefaultBase<Scalar>
1378 IRK1StageTheta_RKBT()
1381 std::ostringstream myDescription;
1382 myDescription << IRK1StageTheta_name() <<
"\n" 1383 <<
"Non-standard finite-difference methods\n" 1384 <<
"in dynamical systems, P. Kama,\n" 1385 <<
"Dissertation, University of Pretoria, pg. 49.\n" 1386 <<
"Comment: Generalized Implicit Midpoint Method\n" 1387 <<
"c = [ theta ]'\n" 1388 <<
"A = [ theta ]\n" 1392 this->setMyDescription(myDescription.str());
1393 typedef ScalarTraits<Scalar> ST;
1394 theta_default_ = ST::one()/(2*ST::one());
1395 theta_ = theta_default_;
1398 RCP<ParameterList> validPL = Teuchos::parameterList();
1399 validPL->set(
"Description",
"",this->getMyDescription());
1400 validPL->set<
double>(
"theta",theta_default_,
1401 "Valid values are 0 <= theta <= 1, where theta = 0 " 1402 "implies Forward Euler, theta = 1/2 implies midpoint " 1403 "method, and theta = 1 implies Backward Euler. " 1404 "For theta != 1/2, this method is first-order accurate, " 1405 "and with theta = 1/2, it is second-order accurate. " 1406 "This method is A-stable, but becomes L-stable with theta=1.");
1407 Teuchos::setupVerboseObjectSublist(&*validPL);
1408 this->setMyValidParameterList(validPL);
1413 typedef ScalarTraits<Scalar> ST;
1414 int myNumStages = 1;
1415 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1416 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1417 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1424 if (theta_ == theta_default_) this->setMy_order(2);
1425 else this->setMy_order(1);
1428 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
1430 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1431 paramList->validateParameters(*this->getValidParameters());
1432 Teuchos::readVerboseObjectSublist(&*paramList,
this);
1433 theta_ = paramList->get<
double>(
"theta",theta_default_);
1435 this->setMyParamList(paramList);
1438 Scalar theta_default_;
1443 template<
class Scalar>
1444 class IRK2StageTheta_RKBT :
1445 virtual public RKButcherTableauDefaultBase<Scalar>
1448 IRK2StageTheta_RKBT()
1450 std::ostringstream myDescription;
1451 myDescription << IRK2StageTheta_name() <<
"\n" 1452 <<
"Computer Methods for ODEs and DAEs\n" 1453 <<
"U. M. Ascher and L. R. Petzold\n" 1457 <<
" [ 1-theta theta ]\n" 1458 <<
"b = [ 1-theta theta ]'\n" 1461 this->setMyDescription(myDescription.str());
1462 typedef ScalarTraits<Scalar> ST;
1463 theta_default_ = ST::one()/(2*ST::one());
1464 theta_ = theta_default_;
1467 RCP<ParameterList> validPL = Teuchos::parameterList();
1468 validPL->set(
"Description",
"",this->getMyDescription());
1469 validPL->set<
double>(
"theta",theta_default_,
1470 "Valid values are 0 < theta <= 1, where theta = 0 " 1471 "implies Forward Euler, theta = 1/2 implies trapezoidal " 1472 "method, and theta = 1 implies Backward Euler. " 1473 "For theta != 1/2, this method is first-order accurate, " 1474 "and with theta = 1/2, it is second-order accurate. " 1475 "This method is A-stable, but becomes L-stable with theta=1.");
1476 Teuchos::setupVerboseObjectSublist(&*validPL);
1477 this->setMyValidParameterList(validPL);
1481 typedef ScalarTraits<Scalar> ST;
1482 int myNumStages = 2;
1483 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1484 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1485 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1486 Scalar one = ST::one();
1487 Scalar zero = ST::zero();
1490 myA(1,0) = as<Scalar>( one - theta_ );
1492 myb(0) = as<Scalar>( one - theta_ );
1500 TEUCHOS_TEST_FOR_EXCEPTION(
1501 theta_ == zero, std::logic_error,
1502 "'theta' can not be zero, as it makes this IRK stepper explicit.");
1503 if (theta_ == theta_default_) this->setMy_order(2);
1504 else this->setMy_order(1);
1506 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
1508 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1509 paramList->validateParameters(*this->getValidParameters());
1510 Teuchos::readVerboseObjectSublist(&*paramList,
this);
1511 theta_ = paramList->get<
double>(
"theta",theta_default_);
1513 this->setMyParamList(paramList);
1516 Scalar theta_default_;
1521 template<
class Scalar>
1522 class Implicit1Stage2ndOrderGauss_RKBT :
1523 virtual public RKButcherTableauDefaultBase<Scalar>
1526 Implicit1Stage2ndOrderGauss_RKBT()
1529 std::ostringstream myDescription;
1530 myDescription << Implicit1Stage2ndOrderGauss_name() <<
"\n" 1532 <<
"Solving Ordinary Differential Equations II:\n" 1533 <<
"Stiff and Differential-Algebraic Problems,\n" 1534 <<
"2nd Revised Edition\n" 1535 <<
"E. Hairer and G. Wanner\n" 1536 <<
"Table 5.2, pg 72\n" 1537 <<
"Also: Implicit midpoint rule\n" 1538 <<
"Solving Ordinary Differential Equations I:\n" 1539 <<
"Nonstiff Problems, 2nd Revised Edition\n" 1540 <<
"E. Hairer, S. P. Norsett, and G. Wanner\n" 1541 <<
"Table 7.1, pg 205\n" 1544 <<
"b = [ 1 ]'" << std::endl;
1545 typedef ScalarTraits<Scalar> ST;
1546 int myNumStages = 1;
1547 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1548 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1549 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1550 Scalar onehalf = ST::one()/(2*ST::one());
1551 Scalar one = ST::one();
1555 this->setMyDescription(myDescription.str());
1559 this->setMy_order(2);
1564 template<
class Scalar>
1565 class Implicit2Stage4thOrderGauss_RKBT :
1566 virtual public RKButcherTableauDefaultBase<Scalar>
1569 Implicit2Stage4thOrderGauss_RKBT()
1572 std::ostringstream myDescription;
1573 myDescription << Implicit2Stage4thOrderGauss_name() <<
"\n" 1575 <<
"Solving Ordinary Differential Equations II:\n" 1576 <<
"Stiff and Differential-Algebraic Problems,\n" 1577 <<
"2nd Revised Edition\n" 1578 <<
"E. Hairer and G. Wanner\n" 1579 <<
"Table 5.2, pg 72\n" 1580 <<
"c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n" 1581 <<
"A = [ 1/4 1/4-sqrt(3)/6 ]\n" 1582 <<
" [ 1/4+sqrt(3)/6 1/4 ]\n" 1583 <<
"b = [ 1/2 1/2 ]'" << std::endl;
1584 typedef ScalarTraits<Scalar> ST;
1585 int myNumStages = 2;
1586 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1587 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1588 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1589 Scalar one = ST::one();
1590 Scalar onehalf = as<Scalar>(one/(2*one));
1591 Scalar three = as<Scalar>(3*one);
1592 Scalar six = as<Scalar>(6*one);
1593 Scalar onefourth = as<Scalar>(one/(4*one));
1594 Scalar alpha = ST::squareroot(three)/six;
1596 myA(0,0) = onefourth;
1597 myA(0,1) = onefourth-alpha;
1598 myA(1,0) = onefourth+alpha;
1599 myA(1,1) = onefourth;
1602 myc(0) = onehalf-alpha;
1603 myc(1) = onehalf+alpha;
1604 this->setMyDescription(myDescription.str());
1608 this->setMy_order(4);
1613 template<
class Scalar>
1614 class Implicit3Stage6thOrderGauss_RKBT :
1615 virtual public RKButcherTableauDefaultBase<Scalar>
1618 Implicit3Stage6thOrderGauss_RKBT()
1621 std::ostringstream myDescription;
1622 myDescription << Implicit3Stage6thOrderGauss_name() <<
"\n" 1624 <<
"Solving Ordinary Differential Equations II:\n" 1625 <<
"Stiff and Differential-Algebraic Problems,\n" 1626 <<
"2nd Revised Edition\n" 1627 <<
"E. Hairer and G. Wanner\n" 1628 <<
"Table 5.2, pg 72\n" 1629 <<
"c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n" 1630 <<
"A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n" 1631 <<
" [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n" 1632 <<
" [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n" 1633 <<
"b = [ 5/18 4/9 5/18 ]'" << std::endl;
1634 typedef ScalarTraits<Scalar> ST;
1635 int myNumStages = 3;
1636 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1637 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1638 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1639 Scalar one = ST::one();
1640 Scalar ten = as<Scalar>(10*one);
1641 Scalar fifteen = as<Scalar>(15*one);
1642 Scalar twentyfour = as<Scalar>(24*one);
1643 Scalar thirty = as<Scalar>(30*one);
1644 Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
1645 Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
1646 Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
1647 Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
1649 myA(0,0) = as<Scalar>(5*one/(36*one));
1650 myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
1651 myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
1652 myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
1653 myA(1,1) = as<Scalar>(2*one/(9*one));
1654 myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
1655 myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
1656 myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
1657 myA(2,2) = as<Scalar>(5*one/(36*one));
1658 myb(0) = as<Scalar>(5*one/(18*one));
1659 myb(1) = as<Scalar>(4*one/(9*one));
1660 myb(2) = as<Scalar>(5*one/(18*one));
1661 myc(0) = as<Scalar>(one/(2*one))-sqrt15over10;
1662 myc(1) = as<Scalar>(one/(2*one));
1663 myc(2) = as<Scalar>(one/(2*one))+sqrt15over10;
1664 this->setMyDescription(myDescription.str());
1668 this->setMy_order(6);
1673 template<
class Scalar>
1674 class Implicit1Stage1stOrderRadauA_RKBT :
1675 virtual public RKButcherTableauDefaultBase<Scalar>
1678 Implicit1Stage1stOrderRadauA_RKBT()
1681 std::ostringstream myDescription;
1682 myDescription << Implicit1Stage1stOrderRadauA_name() <<
"\n" 1684 <<
"Solving Ordinary Differential Equations II:\n" 1685 <<
"Stiff and Differential-Algebraic Problems,\n" 1686 <<
"2nd Revised Edition\n" 1687 <<
"E. Hairer and G. Wanner\n" 1688 <<
"Table 5.3, pg 73\n" 1691 <<
"b = [ 1 ]'" << std::endl;
1692 typedef ScalarTraits<Scalar> ST;
1693 int myNumStages = 1;
1694 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1695 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1696 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1697 Scalar one = ST::one();
1698 Scalar zero = ST::zero();
1702 this->setMyDescription(myDescription.str());
1706 this->setMy_order(1);
1711 template<
class Scalar>
1712 class Implicit2Stage3rdOrderRadauA_RKBT :
1713 virtual public RKButcherTableauDefaultBase<Scalar>
1716 Implicit2Stage3rdOrderRadauA_RKBT()
1719 std::ostringstream myDescription;
1720 myDescription << Implicit2Stage3rdOrderRadauA_name() <<
"\n" 1722 <<
"Solving Ordinary Differential Equations II:\n" 1723 <<
"Stiff and Differential-Algebraic Problems,\n" 1724 <<
"2nd Revised Edition\n" 1725 <<
"E. Hairer and G. Wanner\n" 1726 <<
"Table 5.3, pg 73\n" 1727 <<
"c = [ 0 2/3 ]'\n" 1728 <<
"A = [ 1/4 -1/4 ]\n" 1729 <<
" [ 1/4 5/12 ]\n" 1730 <<
"b = [ 1/4 3/4 ]'" << std::endl;
1731 typedef ScalarTraits<Scalar> ST;
1732 int myNumStages = 2;
1733 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1734 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1735 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1736 Scalar zero = ST::zero();
1737 Scalar one = ST::one();
1738 myA(0,0) = as<Scalar>(one/(4*one));
1739 myA(0,1) = as<Scalar>(-one/(4*one));
1740 myA(1,0) = as<Scalar>(one/(4*one));
1741 myA(1,1) = as<Scalar>(5*one/(12*one));
1742 myb(0) = as<Scalar>(one/(4*one));
1743 myb(1) = as<Scalar>(3*one/(4*one));
1745 myc(1) = as<Scalar>(2*one/(3*one));
1746 this->setMyDescription(myDescription.str());
1750 this->setMy_order(3);
1755 template<
class Scalar>
1756 class Implicit3Stage5thOrderRadauA_RKBT :
1757 virtual public RKButcherTableauDefaultBase<Scalar>
1760 Implicit3Stage5thOrderRadauA_RKBT()
1763 std::ostringstream myDescription;
1764 myDescription << Implicit3Stage5thOrderRadauA_name() <<
"\n" 1766 <<
"Solving Ordinary Differential Equations II:\n" 1767 <<
"Stiff and Differential-Algebraic Problems,\n" 1768 <<
"2nd Revised Edition\n" 1769 <<
"E. Hairer and G. Wanner\n" 1770 <<
"Table 5.4, pg 73\n" 1771 <<
"c = [ 0 (6-sqrt(6))/10 (6+sqrt(6))/10 ]'\n" 1772 <<
"A = [ 1/9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 ]\n" 1773 <<
" [ 1/9 (88+7*sqrt(6))/360 (88-43*sqrt(6))/360 ]\n" 1774 <<
" [ 1/9 (88+43*sqrt(6))/360 (88-7*sqrt(6))/360 ]\n" 1775 <<
"b = [ 1/9 (16+sqrt(6))/36 (16-sqrt(6))/36 ]'" << std::endl;
1776 typedef ScalarTraits<Scalar> ST;
1777 int myNumStages = 3;
1778 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1779 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1780 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1781 Scalar zero = ST::zero();
1782 Scalar one = ST::one();
1783 myA(0,0) = as<Scalar>(one/(9*one));
1784 myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
1785 myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
1786 myA(1,0) = as<Scalar>(one/(9*one));
1787 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1788 myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
1789 myA(2,0) = as<Scalar>(one/(9*one));
1790 myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
1791 myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1792 myb(0) = as<Scalar>(one/(9*one));
1793 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1794 myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1796 myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
1797 myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
1798 this->setMyDescription(myDescription.str());
1802 this->setMy_order(5);
1807 template<
class Scalar>
1808 class Implicit1Stage1stOrderRadauB_RKBT :
1809 virtual public RKButcherTableauDefaultBase<Scalar>
1812 Implicit1Stage1stOrderRadauB_RKBT()
1815 std::ostringstream myDescription;
1816 myDescription << Implicit1Stage1stOrderRadauB_name() <<
"\n" 1818 <<
"Solving Ordinary Differential Equations II:\n" 1819 <<
"Stiff and Differential-Algebraic Problems,\n" 1820 <<
"2nd Revised Edition\n" 1821 <<
"E. Hairer and G. Wanner\n" 1822 <<
"Table 5.5, pg 74\n" 1825 <<
"b = [ 1 ]'" << std::endl;
1826 typedef ScalarTraits<Scalar> ST;
1827 int myNumStages = 1;
1828 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1829 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1830 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1831 Scalar one = ST::one();
1835 this->setMyDescription(myDescription.str());
1839 this->setMy_order(1);
1844 template<
class Scalar>
1845 class Implicit2Stage3rdOrderRadauB_RKBT :
1846 virtual public RKButcherTableauDefaultBase<Scalar>
1849 Implicit2Stage3rdOrderRadauB_RKBT()
1852 std::ostringstream myDescription;
1853 myDescription << Implicit2Stage3rdOrderRadauB_name() <<
"\n" 1855 <<
"Solving Ordinary Differential Equations II:\n" 1856 <<
"Stiff and Differential-Algebraic Problems,\n" 1857 <<
"2nd Revised Edition\n" 1858 <<
"E. Hairer and G. Wanner\n" 1859 <<
"Table 5.5, pg 74\n" 1860 <<
"c = [ 1/3 1 ]'\n" 1861 <<
"A = [ 5/12 -1/12 ]\n" 1863 <<
"b = [ 3/4 1/4 ]'" << std::endl;
1864 typedef ScalarTraits<Scalar> ST;
1865 int myNumStages = 2;
1866 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1867 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1868 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1869 Scalar one = ST::one();
1870 myA(0,0) = as<Scalar>( 5*one/(12*one) );
1871 myA(0,1) = as<Scalar>( -one/(12*one) );
1872 myA(1,0) = as<Scalar>( 3*one/(4*one) );
1873 myA(1,1) = as<Scalar>( one/(4*one) );
1874 myb(0) = as<Scalar>( 3*one/(4*one) );
1875 myb(1) = as<Scalar>( one/(4*one) );
1876 myc(0) = as<Scalar>( one/(3*one) );
1878 this->setMyDescription(myDescription.str());
1882 this->setMy_order(3);
1887 template<
class Scalar>
1888 class Implicit3Stage5thOrderRadauB_RKBT :
1889 virtual public RKButcherTableauDefaultBase<Scalar>
1892 Implicit3Stage5thOrderRadauB_RKBT()
1895 std::ostringstream myDescription;
1896 myDescription << Implicit3Stage5thOrderRadauB_name() <<
"\n" 1898 <<
"Solving Ordinary Differential Equations II:\n" 1899 <<
"Stiff and Differential-Algebraic Problems,\n" 1900 <<
"2nd Revised Edition\n" 1901 <<
"E. Hairer and G. Wanner\n" 1902 <<
"Table 5.6, pg 74\n" 1903 <<
"c = [ (4-sqrt(6))/10 (4+sqrt(6))/10 1 ]'\n" 1904 <<
"A = [ A1 A2 A3 ]\n" 1905 <<
" A1 = [ (88-7*sqrt(6))/360 ]\n" 1906 <<
" [ (296+169*sqrt(6))/1800 ]\n" 1907 <<
" [ (16-sqrt(6))/36 ]\n" 1908 <<
" A2 = [ (296-169*sqrt(6))/1800 ]\n" 1909 <<
" [ (88+7*sqrt(6))/360 ]\n" 1910 <<
" [ (16+sqrt(6))/36 ]\n" 1911 <<
" A3 = [ (-2+3*sqrt(6))/225 ]\n" 1912 <<
" [ (-2-3*sqrt(6))/225 ]\n" 1914 <<
"b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'" 1916 typedef ScalarTraits<Scalar> ST;
1917 int myNumStages = 3;
1918 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1919 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1920 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1921 Scalar one = ST::one();
1922 myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1923 myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
1924 myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
1925 myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
1926 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1927 myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
1928 myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1929 myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1930 myA(2,2) = as<Scalar>( one/(9*one) );
1931 myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1932 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1933 myb(2) = as<Scalar>( one/(9*one) );
1934 myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
1935 myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
1937 this->setMyDescription(myDescription.str());
1941 this->setMy_order(5);
1946 template<
class Scalar>
1947 class Implicit2Stage2ndOrderLobattoA_RKBT :
1948 virtual public RKButcherTableauDefaultBase<Scalar>
1951 Implicit2Stage2ndOrderLobattoA_RKBT()
1954 std::ostringstream myDescription;
1955 myDescription << Implicit2Stage2ndOrderLobattoA_name() <<
"\n" 1957 <<
"Solving Ordinary Differential Equations II:\n" 1958 <<
"Stiff and Differential-Algebraic Problems,\n" 1959 <<
"2nd Revised Edition\n" 1960 <<
"E. Hairer and G. Wanner\n" 1961 <<
"Table 5.7, pg 75\n" 1965 <<
"b = [ 1/2 1/2 ]'" << std::endl;
1966 typedef ScalarTraits<Scalar> ST;
1967 int myNumStages = 2;
1968 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1969 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1970 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1971 Scalar zero = ST::zero();
1972 Scalar one = ST::one();
1975 myA(1,0) = as<Scalar>( one/(2*one) );
1976 myA(1,1) = as<Scalar>( one/(2*one) );
1977 myb(0) = as<Scalar>( one/(2*one) );
1978 myb(1) = as<Scalar>( one/(2*one) );
1981 this->setMyDescription(myDescription.str());
1985 this->setMy_order(2);
1990 template<
class Scalar>
1991 class Implicit3Stage4thOrderLobattoA_RKBT :
1992 virtual public RKButcherTableauDefaultBase<Scalar>
1995 Implicit3Stage4thOrderLobattoA_RKBT()
1998 std::ostringstream myDescription;
1999 myDescription << Implicit3Stage4thOrderLobattoA_name() <<
"\n" 2001 <<
"Solving Ordinary Differential Equations II:\n" 2002 <<
"Stiff and Differential-Algebraic Problems,\n" 2003 <<
"2nd Revised Edition\n" 2004 <<
"E. Hairer and G. Wanner\n" 2005 <<
"Table 5.7, pg 75\n" 2006 <<
"c = [ 0 1/2 1 ]'\n" 2007 <<
"A = [ 0 0 0 ]\n" 2008 <<
" [ 5/24 1/3 -1/24 ]\n" 2009 <<
" [ 1/6 2/3 1/6 ]\n" 2010 <<
"b = [ 1/6 2/3 1/6 ]'" << std::endl;
2011 typedef ScalarTraits<Scalar> ST;
2012 int myNumStages = 3;
2013 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2014 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2015 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2016 Scalar zero = ST::zero();
2017 Scalar one = ST::one();
2021 myA(1,0) = as<Scalar>( (5*one)/(24*one) );
2022 myA(1,1) = as<Scalar>( (one)/(3*one) );
2023 myA(1,2) = as<Scalar>( (-one)/(24*one) );
2024 myA(2,0) = as<Scalar>( (one)/(6*one) );
2025 myA(2,1) = as<Scalar>( (2*one)/(3*one) );
2026 myA(2,2) = as<Scalar>( (1*one)/(6*one) );
2027 myb(0) = as<Scalar>( (one)/(6*one) );
2028 myb(1) = as<Scalar>( (2*one)/(3*one) );
2029 myb(2) = as<Scalar>( (1*one)/(6*one) );
2031 myc(1) = as<Scalar>( one/(2*one) );
2033 this->setMyDescription(myDescription.str());
2037 this->setMy_order(4);
2042 template<
class Scalar>
2043 class Implicit4Stage6thOrderLobattoA_RKBT :
2044 virtual public RKButcherTableauDefaultBase<Scalar>
2047 Implicit4Stage6thOrderLobattoA_RKBT()
2051 typedef Teuchos::ScalarTraits<Scalar> ST;
2053 std::ostringstream myDescription;
2054 myDescription << Implicit4Stage6thOrderLobattoA_name() <<
"\n" 2056 <<
"Solving Ordinary Differential Equations II:\n" 2057 <<
"Stiff and Differential-Algebraic Problems,\n" 2058 <<
"2nd Revised Edition\n" 2059 <<
"E. Hairer and G. Wanner\n" 2060 <<
"Table 5.8, pg 75\n" 2061 <<
"c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n" 2062 <<
"A = [ A1 A2 A3 A4 ]\n" 2064 <<
" [ (11+sqrt(5)/120 ]\n" 2065 <<
" [ (11-sqrt(5)/120 ]\n" 2068 <<
" [ (25-sqrt(5))/120 ]\n" 2069 <<
" [ (25+13*sqrt(5))/120 ]\n" 2072 <<
" [ (25-13*sqrt(5))/120 ]\n" 2073 <<
" [ (25+sqrt(5))/120 ]\n" 2076 <<
" [ (-1+sqrt(5))/120 ]\n" 2077 <<
" [ (-1-sqrt(5))/120 ]\n" 2079 <<
"b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2080 typedef ScalarTraits<Scalar> ST;
2081 int myNumStages = 4;
2082 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2083 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2084 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2085 Scalar zero = ST::zero();
2086 Scalar one = ST::one();
2091 myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
2092 myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
2093 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
2094 myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
2095 myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
2096 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
2097 myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
2098 myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
2099 myA(3,0) = as<Scalar>( one/(12*one) );
2100 myA(3,1) = as<Scalar>( 5*one/(12*one) );
2101 myA(3,2) = as<Scalar>( 5*one/(12*one) );
2102 myA(3,3) = as<Scalar>( one/(12*one) );
2103 myb(0) = as<Scalar>( one/(12*one) );
2104 myb(1) = as<Scalar>( 5*one/(12*one) );
2105 myb(2) = as<Scalar>( 5*one/(12*one) );
2106 myb(3) = as<Scalar>( one/(12*one) );
2108 myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
2109 myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
2111 this->setMyDescription(myDescription.str());
2115 this->setMy_order(6);
2120 template<
class Scalar>
2121 class Implicit2Stage2ndOrderLobattoB_RKBT :
2122 virtual public RKButcherTableauDefaultBase<Scalar>
2125 Implicit2Stage2ndOrderLobattoB_RKBT()
2128 std::ostringstream myDescription;
2129 myDescription << Implicit2Stage2ndOrderLobattoB_name() <<
"\n" 2131 <<
"Solving Ordinary Differential Equations II:\n" 2132 <<
"Stiff and Differential-Algebraic Problems,\n" 2133 <<
"2nd Revised Edition\n" 2134 <<
"E. Hairer and G. Wanner\n" 2135 <<
"Table 5.9, pg 76\n" 2137 <<
"A = [ 1/2 0 ]\n" 2139 <<
"b = [ 1/2 1/2 ]'" << std::endl;
2140 typedef ScalarTraits<Scalar> ST;
2141 int myNumStages = 2;
2142 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2143 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2144 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2145 Scalar zero = ST::zero();
2146 Scalar one = ST::one();
2147 myA(0,0) = as<Scalar>( one/(2*one) );
2149 myA(1,0) = as<Scalar>( one/(2*one) );
2151 myb(0) = as<Scalar>( one/(2*one) );
2152 myb(1) = as<Scalar>( one/(2*one) );
2155 this->setMyDescription(myDescription.str());
2159 this->setMy_order(2);
2164 template<
class Scalar>
2165 class Implicit3Stage4thOrderLobattoB_RKBT :
2166 virtual public RKButcherTableauDefaultBase<Scalar>
2169 Implicit3Stage4thOrderLobattoB_RKBT()
2172 std::ostringstream myDescription;
2173 myDescription << Implicit3Stage4thOrderLobattoB_name() <<
"\n" 2175 <<
"Solving Ordinary Differential Equations II:\n" 2176 <<
"Stiff and Differential-Algebraic Problems,\n" 2177 <<
"2nd Revised Edition\n" 2178 <<
"E. Hairer and G. Wanner\n" 2179 <<
"Table 5.9, pg 76\n" 2180 <<
"c = [ 0 1/2 1 ]'\n" 2181 <<
"A = [ 1/6 -1/6 0 ]\n" 2182 <<
" [ 1/6 1/3 0 ]\n" 2183 <<
" [ 1/6 5/6 0 ]\n" 2184 <<
"b = [ 1/6 2/3 1/6 ]'" << std::endl;
2185 typedef ScalarTraits<Scalar> ST;
2186 int myNumStages = 3;
2187 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2188 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2189 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2190 Scalar zero = ST::zero();
2191 Scalar one = ST::one();
2192 myA(0,0) = as<Scalar>( one/(6*one) );
2193 myA(0,1) = as<Scalar>( -one/(6*one) );
2195 myA(1,0) = as<Scalar>( one/(6*one) );
2196 myA(1,1) = as<Scalar>( one/(3*one) );
2198 myA(2,0) = as<Scalar>( one/(6*one) );
2199 myA(2,1) = as<Scalar>( 5*one/(6*one) );
2201 myb(0) = as<Scalar>( one/(6*one) );
2202 myb(1) = as<Scalar>( 2*one/(3*one) );
2203 myb(2) = as<Scalar>( one/(6*one) );
2205 myc(1) = as<Scalar>( one/(2*one) );
2207 this->setMyDescription(myDescription.str());
2211 this->setMy_order(4);
2216 template<
class Scalar>
2217 class Implicit4Stage6thOrderLobattoB_RKBT :
2218 virtual public RKButcherTableauDefaultBase<Scalar>
2221 Implicit4Stage6thOrderLobattoB_RKBT()
2224 std::ostringstream myDescription;
2225 myDescription << Implicit4Stage6thOrderLobattoB_name() <<
"\n" 2227 <<
"Solving Ordinary Differential Equations II:\n" 2228 <<
"Stiff and Differential-Algebraic Problems,\n" 2229 <<
"2nd Revised Edition\n" 2230 <<
"E. Hairer and G. Wanner\n" 2231 <<
"Table 5.10, pg 76\n" 2232 <<
"c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n" 2233 <<
"A = [ 1/12 (-1-sqrt(5))/24 (-1+sqrt(5))/24 0 ]\n" 2234 <<
" [ 1/12 (25+sqrt(5))/120 (25-13*sqrt(5))/120 0 ]\n" 2235 <<
" [ 1/12 (25+13*sqrt(5))/120 (25-sqrt(5))/120 0 ]\n" 2236 <<
" [ 1/12 (11-sqrt(5))/24 (11+sqrt(5))/24 0 ]\n" 2237 <<
"b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2238 typedef ScalarTraits<Scalar> ST;
2239 int myNumStages = 4;
2240 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2241 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2242 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2243 Scalar zero = ST::zero();
2244 Scalar one = ST::one();
2245 myA(0,0) = as<Scalar>( one/(12*one) );
2246 myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
2247 myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
2249 myA(1,0) = as<Scalar>( one/(12*one) );
2250 myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
2251 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
2253 myA(2,0) = as<Scalar>( one/(12*one) );
2254 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
2255 myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
2257 myA(3,0) = as<Scalar>( one/(12*one) );
2258 myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
2259 myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
2261 myb(0) = as<Scalar>( one/(12*one) );
2262 myb(1) = as<Scalar>( 5*one/(12*one) );
2263 myb(2) = as<Scalar>( 5*one/(12*one) );
2264 myb(3) = as<Scalar>( one/(12*one) );
2266 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2267 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2269 this->setMyDescription(myDescription.str());
2273 this->setMy_order(6);
2278 template<
class Scalar>
2279 class Implicit2Stage2ndOrderLobattoC_RKBT :
2280 virtual public RKButcherTableauDefaultBase<Scalar>
2283 Implicit2Stage2ndOrderLobattoC_RKBT()
2286 std::ostringstream myDescription;
2287 myDescription << Implicit2Stage2ndOrderLobattoC_name() <<
"\n" 2289 <<
"Solving Ordinary Differential Equations II:\n" 2290 <<
"Stiff and Differential-Algebraic Problems,\n" 2291 <<
"2nd Revised Edition\n" 2292 <<
"E. Hairer and G. Wanner\n" 2293 <<
"Table 5.11, pg 76\n" 2295 <<
"A = [ 1/2 -1/2 ]\n" 2297 <<
"b = [ 1/2 1/2 ]'" << std::endl;
2298 typedef ScalarTraits<Scalar> ST;
2299 int myNumStages = 2;
2300 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2301 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2302 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2303 Scalar zero = ST::zero();
2304 Scalar one = ST::one();
2305 myA(0,0) = as<Scalar>( one/(2*one) );
2306 myA(0,1) = as<Scalar>( -one/(2*one) );
2307 myA(1,0) = as<Scalar>( one/(2*one) );
2308 myA(1,1) = as<Scalar>( one/(2*one) );
2309 myb(0) = as<Scalar>( one/(2*one) );
2310 myb(1) = as<Scalar>( one/(2*one) );
2313 this->setMyDescription(myDescription.str());
2317 this->setMy_order(2);
2322 template<
class Scalar>
2323 class Implicit3Stage4thOrderLobattoC_RKBT :
2324 virtual public RKButcherTableauDefaultBase<Scalar>
2327 Implicit3Stage4thOrderLobattoC_RKBT()
2330 std::ostringstream myDescription;
2331 myDescription << Implicit3Stage4thOrderLobattoC_name() <<
"\n" 2333 <<
"Solving Ordinary Differential Equations II:\n" 2334 <<
"Stiff and Differential-Algebraic Problems,\n" 2335 <<
"2nd Revised Edition\n" 2336 <<
"E. Hairer and G. Wanner\n" 2337 <<
"Table 5.11, pg 76\n" 2338 <<
"c = [ 0 1/2 1 ]'\n" 2339 <<
"A = [ 1/6 -1/3 1/6 ]\n" 2340 <<
" [ 1/6 5/12 -1/12 ]\n" 2341 <<
" [ 1/6 2/3 1/6 ]\n" 2342 <<
"b = [ 1/6 2/3 1/6 ]'" << std::endl;
2343 typedef ScalarTraits<Scalar> ST;
2344 int myNumStages = 3;
2345 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2346 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2347 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2348 Scalar zero = ST::zero();
2349 Scalar one = ST::one();
2350 myA(0,0) = as<Scalar>( one/(6*one) );
2351 myA(0,1) = as<Scalar>( -one/(3*one) );
2352 myA(0,2) = as<Scalar>( one/(6*one) );
2353 myA(1,0) = as<Scalar>( one/(6*one) );
2354 myA(1,1) = as<Scalar>( 5*one/(12*one) );
2355 myA(1,2) = as<Scalar>( -one/(12*one) );
2356 myA(2,0) = as<Scalar>( one/(6*one) );
2357 myA(2,1) = as<Scalar>( 2*one/(3*one) );
2358 myA(2,2) = as<Scalar>( one/(6*one) );
2359 myb(0) = as<Scalar>( one/(6*one) );
2360 myb(1) = as<Scalar>( 2*one/(3*one) );
2361 myb(2) = as<Scalar>( one/(6*one) );
2363 myc(1) = as<Scalar>( one/(2*one) );
2365 this->setMyDescription(myDescription.str());
2369 this->setMy_order(4);
2374 template<
class Scalar>
2375 class Implicit4Stage6thOrderLobattoC_RKBT :
2376 virtual public RKButcherTableauDefaultBase<Scalar>
2379 Implicit4Stage6thOrderLobattoC_RKBT()
2382 std::ostringstream myDescription;
2383 myDescription << Implicit4Stage6thOrderLobattoC_name() <<
"\n" 2385 <<
"Solving Ordinary Differential Equations II:\n" 2386 <<
"Stiff and Differential-Algebraic Problems,\n" 2387 <<
"2nd Revised Edition\n" 2388 <<
"E. Hairer and G. Wanner\n" 2389 <<
"Table 5.12, pg 76\n" 2390 <<
"c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n" 2391 <<
"A = [ 1/12 -sqrt(5)/12 sqrt(5)/12 -1/12 ]\n" 2392 <<
" [ 1/12 1/4 (10-7*sqrt(5))/60 sqrt(5)/60 ]\n" 2393 <<
" [ 1/12 (10+7*sqrt(5))/60 1/4 -sqrt(5)/60 ]\n" 2394 <<
" [ 1/12 5/12 5/12 1/12 ]\n" 2395 <<
"b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2396 typedef ScalarTraits<Scalar> ST;
2397 int myNumStages = 4;
2398 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2399 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2400 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2401 Scalar zero = ST::zero();
2402 Scalar one = ST::one();
2403 myA(0,0) = as<Scalar>( one/(12*one) );
2404 myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
2405 myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
2406 myA(0,3) = as<Scalar>( -one/(12*one) );
2407 myA(1,0) = as<Scalar>( one/(12*one) );
2408 myA(1,1) = as<Scalar>( one/(4*one) );
2409 myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
2410 myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
2411 myA(2,0) = as<Scalar>( one/(12*one) );
2412 myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
2413 myA(2,2) = as<Scalar>( one/(4*one) );
2414 myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
2415 myA(3,0) = as<Scalar>( one/(12*one) );
2416 myA(3,1) = as<Scalar>( 5*one/(12*one) );
2417 myA(3,2) = as<Scalar>( 5*one/(12*one) );
2418 myA(3,3) = as<Scalar>( one/(12*one) );
2419 myb(0) = as<Scalar>( one/(12*one) );
2420 myb(1) = as<Scalar>( 5*one/(12*one) );
2421 myb(2) = as<Scalar>( 5*one/(12*one) );
2422 myb(3) = as<Scalar>( one/(12*one) );
2424 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2425 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2427 this->setMyDescription(myDescription.str());
2431 this->setMy_order(6);
2437 template<
class Scalar>
2438 class SDIRK5Stage5thOrder_RKBT :
2439 virtual public RKButcherTableauDefaultBase<Scalar>
2442 SDIRK5Stage5thOrder_RKBT()
2445 std::ostringstream myDescription;
2446 myDescription << SDIRK5Stage5thOrder_name() <<
"\n" 2448 <<
"Solving Ordinary Differential Equations II:\n" 2449 <<
"Stiff and Differential-Algebraic Problems,\n" 2450 <<
"2nd Revised Edition\n" 2451 <<
"E. Hairer and G. Wanner\n" 2453 <<
"c = [ (6-sqrt(6))/10 ]\n" 2454 <<
" [ (6+9*sqrt(6))/35 ]\n" 2456 <<
" [ (4-sqrt(6))/10 ]\n" 2457 <<
" [ (4+sqrt(6))/10 ]\n" 2458 <<
"A = [ A1 A2 A3 A4 A5 ]\n" 2459 <<
" A1 = [ (6-sqrt(6))/10 ]\n" 2460 <<
" [ (-6+5*sqrt(6))/14 ]\n" 2461 <<
" [ (888+607*sqrt(6))/2850 ]\n" 2462 <<
" [ (3153-3082*sqrt(6))/14250 ]\n" 2463 <<
" [ (-32583+14638*sqrt(6))/71250 ]\n" 2465 <<
" [ (6-sqrt(6))/10 ]\n" 2466 <<
" [ (126-161*sqrt(6))/1425 ]\n" 2467 <<
" [ (3213+1148*sqrt(6))/28500 ]\n" 2468 <<
" [ (-17199+364*sqrt(6))/142500 ]\n" 2471 <<
" [ (6-sqrt(6))/10 ]\n" 2472 <<
" [ (-267+88*sqrt(6))/500 ]\n" 2473 <<
" [ (1329-544*sqrt(6))/2500 ]\n" 2477 <<
" [ (6-sqrt(6))/10 ]\n" 2478 <<
" [ (-96+131*sqrt(6))/625 ]\n" 2483 <<
" [ (6-sqrt(6))/10 ]\n" 2487 <<
" [ (16-sqrt(6))/36 ]\n" 2488 <<
" [ (16+sqrt(6))/36 ]" << std::endl;
2489 typedef ScalarTraits<Scalar> ST;
2490 int myNumStages = 5;
2491 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2492 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2493 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2494 Scalar zero = ST::zero();
2495 Scalar one = ST::one();
2496 Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
2497 Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) );
2504 myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
2510 myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
2511 myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
2516 myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
2517 myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
2518 myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
2522 myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
2523 myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
2524 myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
2525 myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
2530 myb(2) = as<Scalar>( one/(9*one) );
2531 myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
2532 myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
2535 myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
2537 myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
2538 myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
2540 this->setMyDescription(myDescription.str());
2544 this->setMy_order(5);
2549 template<
class Scalar>
2550 class SDIRK5Stage4thOrder_RKBT :
2551 virtual public RKButcherTableauDefaultBase<Scalar>
2554 SDIRK5Stage4thOrder_RKBT()
2557 std::ostringstream myDescription;
2558 myDescription << SDIRK5Stage4thOrder_name() <<
"\n" 2560 <<
"Solving Ordinary Differential Equations II:\n" 2561 <<
"Stiff and Differential-Algebraic Problems,\n" 2562 <<
"2nd Revised Edition\n" 2563 <<
"E. Hairer and G. Wanner\n" 2565 <<
"c = [ 1/4 3/4 11/20 1/2 1 ]'\n" 2568 <<
" [ 17/50 -1/25 1/4 ]\n" 2569 <<
" [ 371/1360 -137/2720 15/544 1/4 ]\n" 2570 <<
" [ 25/24 -49/48 125/16 -85/12 1/4 ]\n" 2571 <<
"b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'\n" 2572 <<
"b' = [ 59/48 -17/96 225/32 -85/12 0 ]'" << std::endl;
2573 typedef ScalarTraits<Scalar> ST;
2574 int myNumStages = 5;
2575 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2576 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2577 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2578 Scalar zero = ST::zero();
2579 Scalar one = ST::one();
2580 Scalar onequarter = as<Scalar>( one/(4*one) );
2581 myA(0,0) = onequarter;
2587 myA(1,0) = as<Scalar>( one / (2*one) );
2588 myA(1,1) = onequarter;
2593 myA(2,0) = as<Scalar>( 17*one/(50*one) );
2594 myA(2,1) = as<Scalar>( -one/(25*one) );
2595 myA(2,2) = onequarter;
2599 myA(3,0) = as<Scalar>( 371*one/(1360*one) );
2600 myA(3,1) = as<Scalar>( -137*one/(2720*one) );
2601 myA(3,2) = as<Scalar>( 15*one/(544*one) );
2602 myA(3,3) = onequarter;
2605 myA(4,0) = as<Scalar>( 25*one/(24*one) );
2606 myA(4,1) = as<Scalar>( -49*one/(48*one) );
2607 myA(4,2) = as<Scalar>( 125*one/(16*one) );
2608 myA(4,3) = as<Scalar>( -85*one/(12*one) );
2609 myA(4,4) = onequarter;
2611 myb(0) = as<Scalar>( 25*one/(24*one) );
2612 myb(1) = as<Scalar>( -49*one/(48*one) );
2613 myb(2) = as<Scalar>( 125*one/(16*one) );
2614 myb(3) = as<Scalar>( -85*one/(12*one) );
2615 myb(4) = onequarter;
2625 myc(0) = onequarter;
2626 myc(1) = as<Scalar>( 3*one/(4*one) );
2627 myc(2) = as<Scalar>( 11*one/(20*one) );
2628 myc(3) = as<Scalar>( one/(2*one) );
2631 this->setMyDescription(myDescription.str());
2635 this->setMy_order(4);
2640 template<
class Scalar>
2641 class SDIRK3Stage4thOrder_RKBT :
2642 virtual public RKButcherTableauDefaultBase<Scalar>
2645 SDIRK3Stage4thOrder_RKBT()
2648 std::ostringstream myDescription;
2649 myDescription << SDIRK3Stage4thOrder_name() <<
"\n" 2651 <<
"Solving Ordinary Differential Equations II:\n" 2652 <<
"Stiff and Differential-Algebraic Problems,\n" 2653 <<
"2nd Revised Edition\n" 2654 <<
"E. Hairer and G. Wanner\n" 2656 <<
"gamma = (1/sqrt(3))*cos(pi/18)+1/2\n" 2657 <<
"delta = 1/(6*(2*gamma-1)^2)\n" 2658 <<
"c = [ gamma 1/2 1-gamma ]'\n" 2659 <<
"A = [ gamma ]\n" 2660 <<
" [ 1/2-gamma gamma ]\n" 2661 <<
" [ 2*gamma 1-4*gamma gamma ]\n" 2662 <<
"b = [ delta 1-2*delta delta ]'" << std::endl;
2663 typedef ScalarTraits<Scalar> ST;
2664 int myNumStages = 3;
2665 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2666 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2667 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2668 Scalar zero = ST::zero();
2669 Scalar one = ST::one();
2670 Scalar pi = as<Scalar>(4*one)*std::atan(one);
2671 Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
2672 Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
2677 myA(1,0) = as<Scalar>( one/(2*one) - gamma );
2681 myA(2,0) = as<Scalar>( 2*gamma );
2682 myA(2,1) = as<Scalar>( one - 4*gamma );
2686 myb(1) = as<Scalar>( one-2*delta );
2690 myc(1) = as<Scalar>( one/(2*one) );
2691 myc(2) = as<Scalar>( one - gamma );
2693 this->setMyDescription(myDescription.str());
2697 this->setMy_order(4);
2705 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP