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 Explicit3Stage3rdOrderHeun_name() {
return "Explicit 3 Stage 3rd order by Heun"; }
64 inline const std::string Explicit3Stage3rdOrder_name() {
return "Explicit 3 Stage 3rd order"; }
65 inline const std::string Explicit3Stage3rdOrderTVD_name() {
return "Explicit 3 Stage 3rd order TVD"; }
66 inline const std::string Explicit4Stage3rdOrderRunge_name() {
return "Explicit 4 Stage 3rd order by Runge"; }
68 inline const std::string IRK1StageTheta_name() {
return "IRK 1 Stage Theta Method"; }
69 inline const std::string IRK2StageTheta_name() {
return "IRK 2 Stage Theta Method"; }
70 inline const std::string Implicit1Stage2ndOrderGauss_name() {
return "Implicit 1 Stage 2nd order Gauss"; }
71 inline const std::string Implicit2Stage4thOrderGauss_name() {
return "Implicit 2 Stage 4th order Gauss"; }
72 inline const std::string Implicit3Stage6thOrderGauss_name() {
return "Implicit 3 Stage 6th order Gauss"; }
74 inline const std::string Implicit1Stage1stOrderRadauA_name() {
return "Implicit 1 Stage 1st order Radau left"; }
75 inline const std::string Implicit2Stage3rdOrderRadauA_name() {
return "Implicit 2 Stage 3rd order Radau left"; }
76 inline const std::string Implicit3Stage5thOrderRadauA_name() {
return "Implicit 3 Stage 5th order Radau left"; }
78 inline const std::string Implicit1Stage1stOrderRadauB_name() {
return "Implicit 1 Stage 1st order Radau right"; }
79 inline const std::string Implicit2Stage3rdOrderRadauB_name() {
return "Implicit 2 Stage 3rd order Radau right"; }
80 inline const std::string Implicit3Stage5thOrderRadauB_name() {
return "Implicit 3 Stage 5th order Radau right"; }
82 inline const std::string Implicit2Stage2ndOrderLobattoA_name() {
return "Implicit 2 Stage 2nd order Lobatto A"; }
83 inline const std::string Implicit3Stage4thOrderLobattoA_name() {
return "Implicit 3 Stage 4th order Lobatto A"; }
84 inline const std::string Implicit4Stage6thOrderLobattoA_name() {
return "Implicit 4 Stage 6th order Lobatto A"; }
86 inline const std::string Implicit2Stage2ndOrderLobattoB_name() {
return "Implicit 2 Stage 2nd order Lobatto B"; }
87 inline const std::string Implicit3Stage4thOrderLobattoB_name() {
return "Implicit 3 Stage 4th order Lobatto B"; }
88 inline const std::string Implicit4Stage6thOrderLobattoB_name() {
return "Implicit 4 Stage 6th order Lobatto B"; }
90 inline const std::string Implicit2Stage2ndOrderLobattoC_name() {
return "Implicit 2 Stage 2nd order Lobatto C"; }
91 inline const std::string Implicit3Stage4thOrderLobattoC_name() {
return "Implicit 3 Stage 4th order Lobatto C"; }
92 inline const std::string Implicit4Stage6thOrderLobattoC_name() {
return "Implicit 4 Stage 6th order Lobatto C"; }
94 inline const std::string Implicit2Stage4thOrderHammerHollingsworth_name() {
return "Implicit 2 Stage 4th Order Hammer & Hollingsworth"; }
95 inline const std::string Implicit3Stage6thOrderKuntzmannButcher_name() {
return "Implicit 3 Stage 6th Order Kuntzmann & Butcher"; }
96 inline const std::string Implicit4Stage8thOrderKuntzmannButcher_name() {
return "Implicit 4 Stage 8th Order Kuntzmann & Butcher"; }
98 inline const std::string DIRK2Stage3rdOrder_name() {
return "Diagonal IRK 2 Stage 3rd order"; }
100 inline const std::string SDIRK2Stage2ndOrder_name() {
return "Singly Diagonal IRK 2 Stage 2nd order"; }
101 inline const std::string SDIRK2Stage3rdOrder_name() {
return "Singly Diagonal IRK 2 Stage 3rd order"; }
102 inline const std::string SDIRK5Stage5thOrder_name() {
return "Singly Diagonal IRK 5 Stage 5th order"; }
103 inline const std::string SDIRK5Stage4thOrder_name() {
return "Singly Diagonal IRK 5 Stage 4th order"; }
104 inline const std::string SDIRK3Stage4thOrder_name() {
return "Singly Diagonal IRK 3 Stage 4th order"; }
106 template<
class Scalar>
107 class RKButcherTableauDefaultBase :
108 virtual public RKButcherTableauBase<Scalar>,
109 virtual public Teuchos::ParameterListAcceptorDefaultBase
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; }
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,
131 const std::string& longDescription_in
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 );
144 longDescription_ = longDescription_in;
151 virtual void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
153 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
154 paramList->validateParameters(*this->getValidParameters());
155 Teuchos::readVerboseObjectSublist(&*paramList,
this);
156 setMyParamList(paramList);
160 virtual RCP<const Teuchos::ParameterList> getValidParameters()
const 162 if (is_null(validPL_)) {
163 validPL_ = Teuchos::parameterList();
164 validPL_->set(
"Description",
"",this->getMyDescription());
165 Teuchos::setupVerboseObjectSublist(&*validPL_);
176 virtual std::string description()
const {
return "Rythmos::RKButcherTableauDefaultBase"; }
179 virtual void describe(
180 Teuchos::FancyOStream &out,
181 const Teuchos::EVerbosityLevel verbLevel
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;
198 void setMyDescription(std::string longDescription) { longDescription_ = longDescription; }
199 const std::string& getMyDescription()
const {
return longDescription_; }
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; }
206 void setMyValidParameterList(
const RCP<ParameterList> validPL ) { validPL_ = validPL; }
207 RCP<ParameterList> getMyNonconstValidParameterList() {
return validPL_; }
210 Teuchos::SerialDenseMatrix<int,Scalar> A_;
211 Teuchos::SerialDenseVector<int,Scalar> b_;
212 Teuchos::SerialDenseVector<int,Scalar> c_;
214 std::string longDescription_;
215 mutable RCP<ParameterList> validPL_;
220 template<
class Scalar>
221 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
223 return(rcp(
new RKButcherTableauDefaultBase<Scalar>()));
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,
233 const std::string& description_in =
"" 236 RCP<RKButcherTableauDefaultBase<Scalar> > rkbt = rcp(
new RKButcherTableauDefaultBase<Scalar>());
237 rkbt->initialize(A_in,b_in,c_in,order_in,description_in);
242 template<
class Scalar>
243 class BackwardEuler_RKBT :
244 virtual public RKButcherTableauDefaultBase<Scalar>
249 std::ostringstream myDescription;
250 myDescription << RKBT_BackwardEuler_name() <<
"\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);
259 Teuchos::SerialDenseVector<int,Scalar> myc(1);
262 this->setMyDescription(myDescription.str());
266 this->setMy_order(1);
272 template<
class Scalar>
273 class ForwardEuler_RKBT :
274 virtual public RKButcherTableauDefaultBase<Scalar>
280 std::ostringstream myDescription;
281 myDescription << RKBT_ForwardEuler_name() <<
"\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);
289 Teuchos::SerialDenseVector<int,Scalar> myc(1);
291 this->setMyDescription(myDescription.str());
295 this->setMy_order(1);
300 template<
class Scalar>
301 class Explicit4Stage4thOrder_RKBT :
302 virtual public RKButcherTableauDefaultBase<Scalar>
305 Explicit4Stage4thOrder_RKBT()
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" 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());
328 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
329 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
330 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
365 this->setMyDescription(myDescription.str());
369 this->setMy_order(4);
374 template<
class Scalar>
375 class Explicit3_8Rule_RKBT :
376 virtual public RKButcherTableauDefaultBase<Scalar>
379 Explicit3_8Rule_RKBT()
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" 393 <<
"b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
394 typedef ScalarTraits<Scalar> ST;
396 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
397 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
398 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
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()));
413 myA(1,0) = one_third;
418 myA(2,0) = as<Scalar>(-one_third);
424 myA(3,1) = as<Scalar>(-one);
430 myb(1) = three_eighth;
431 myb(2) = three_eighth;
440 this->setMyDescription(myDescription.str());
444 this->setMy_order(4);
449 template<
class Scalar>
450 class Explicit4Stage3rdOrderRunge_RKBT :
451 virtual public RKButcherTableauDefaultBase<Scalar>
454 Explicit4Stage3rdOrderRunge_RKBT()
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" 468 <<
"b = [ 1/6 2/3 0 1/6 ]'" << std::endl;
469 typedef ScalarTraits<Scalar> ST;
471 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
472 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
473 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
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();
514 this->setMyDescription(myDescription.str());
518 this->setMy_order(3);
523 template<
class Scalar>
524 class Explicit3Stage3rdOrder_RKBT :
525 virtual public RKButcherTableauDefaultBase<Scalar>
528 Explicit3Stage3rdOrder_RKBT()
531 std::ostringstream myDescription;
532 myDescription << Explicit3Stage3rdOrder_name() <<
"\n" 533 <<
"c = [ 0 1/2 1 ]'\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());
547 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
548 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
549 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
574 this->setMyDescription(myDescription.str());
578 this->setMy_order(3);
583 template<
class Scalar>
584 class Explicit3Stage3rdOrderTVD_RKBT :
585 virtual public RKButcherTableauDefaultBase<Scalar>
588 Explicit3Stage3rdOrderTVD_RKBT()
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" 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" 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());
616 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
617 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
618 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
629 myA(2,0) = onefourth;
630 myA(2,1) = onefourth;
643 this->setMyDescription(myDescription.str());
647 this->setMy_order(3);
652 template<
class Scalar>
653 class Explicit3Stage3rdOrderHeun_RKBT :
654 virtual public RKButcherTableauDefaultBase<Scalar>
657 Explicit3Stage3rdOrderHeun_RKBT()
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" 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);
679 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
680 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
681 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
693 myA(2,1) = twothirds;
699 myb(2) = threefourths;
706 this->setMyDescription(myDescription.str());
710 this->setMy_order(3);
715 template<
class Scalar>
716 class Explicit2Stage2ndOrderRunge_RKBT :
717 virtual public RKButcherTableauDefaultBase<Scalar>
720 Explicit2Stage2ndOrderRunge_RKBT()
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" 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());
739 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
740 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
741 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
758 this->setMyDescription(myDescription.str());
762 this->setMy_order(2);
767 template<
class Scalar>
768 class ExplicitTrapezoidal_RKBT :
769 virtual public RKButcherTableauDefaultBase<Scalar>
772 ExplicitTrapezoidal_RKBT()
774 std::ostringstream myDescription;
775 myDescription << ExplicitTrapezoidal_name() <<
"\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());
786 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
787 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
788 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
805 this->setMyDescription(myDescription.str());
809 this->setMy_order(2);
814 template<
class Scalar>
815 class SDIRK2Stage2ndOrder_RKBT :
816 virtual public RKButcherTableauDefaultBase<Scalar>
819 SDIRK2Stage2ndOrder_RKBT()
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" 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;
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_;
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);
851 typedef ScalarTraits<Scalar> ST;
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();
860 myA(1,0) = as<Scalar>( one - gamma_ );
862 myb(0) = as<Scalar>( one - gamma_ );
870 this->setMy_order(2);
872 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
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_);
879 this->setMyParamList(paramList);
882 Scalar gamma_default_;
889 template<
class Scalar>
890 class SDIRK2Stage3rdOrder_RKBT :
891 virtual public RKButcherTableauDefaultBase<Scalar>
894 SDIRK2Stage3rdOrder_RKBT()
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;
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_;
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.");
935 Teuchos::setupVerboseObjectSublist(&*validPL);
936 this->setMyValidParameterList(validPL);
940 typedef ScalarTraits<Scalar> ST;
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();
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) );
956 myA(1,0) = as<Scalar>( one - 2*gamma );
958 myb(0) = as<Scalar>( one/(2*one) );
959 myb(1) = as<Scalar>( one/(2*one) );
961 myc(1) = as<Scalar>( one - gamma );
966 this->setMy_order(3);
968 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
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_);
983 this->setMyParamList(paramList);
986 bool thirdOrderAStable_default_;
987 bool thirdOrderAStable_;
988 bool secondOrderLStable_default_;
989 bool secondOrderLStable_;
990 Scalar gamma_default_;
995 template<
class Scalar>
996 class DIRK2Stage3rdOrder_RKBT :
997 virtual public RKButcherTableauDefaultBase<Scalar>
1000 DIRK2Stage3rdOrder_RKBT()
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" 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();
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) );
1028 myc(1) = as<Scalar>( 2*one/(3*one) );
1029 this->setMyDescription(myDescription.str());
1033 this->setMy_order(3);
1038 template<
class Scalar>
1039 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
1040 virtual public RKButcherTableauDefaultBase<Scalar>
1043 Implicit3Stage6thOrderKuntzmannButcher_RKBT()
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());
1083 this->setMy_order(6);
1088 template<
class Scalar>
1089 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
1090 virtual public RKButcherTableauDefaultBase<Scalar>
1093 Implicit4Stage8thOrderKuntzmannButcher_RKBT()
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" 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 );
1137 myA(0,1) = w1p-w3+w4p;
1138 myA(0,2) = w1p-w3-w4p;
1140 myA(1,0) = w1-w3p+w4;
1143 myA(1,3) = w1-w3p-w4;
1144 myA(2,0) = w1+w3p+w4;
1147 myA(2,3) = w1+w3p-w4;
1149 myA(3,1) = w1p+w3+w4p;
1150 myA(3,2) = w1p+w3-w4p;
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());
1164 this->setMy_order(8);
1169 template<
class Scalar>
1170 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
1171 virtual public RKButcherTableauDefaultBase<Scalar>
1174 Implicit2Stage4thOrderHammerHollingsworth_RKBT()
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;
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());
1208 this->setMy_order(4);
1213 template<
class Scalar>
1214 class IRK1StageTheta_RKBT :
1215 virtual public RKButcherTableauDefaultBase<Scalar>
1218 IRK1StageTheta_RKBT()
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" 1232 this->setMyDescription(myDescription.str());
1233 typedef ScalarTraits<Scalar> ST;
1234 theta_default_ = ST::one()/(2*ST::one());
1235 theta_ = theta_default_;
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);
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);
1264 if (theta_ == theta_default_) this->setMy_order(2);
1265 else this->setMy_order(1);
1268 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
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_);
1275 this->setMyParamList(paramList);
1278 Scalar theta_default_;
1283 template<
class Scalar>
1284 class IRK2StageTheta_RKBT :
1285 virtual public RKButcherTableauDefaultBase<Scalar>
1288 IRK2StageTheta_RKBT()
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" 1297 <<
" [ 1-theta theta ]\n" 1298 <<
"b = [ 1-theta theta ]'\n" 1301 this->setMyDescription(myDescription.str());
1302 typedef ScalarTraits<Scalar> ST;
1303 theta_default_ = ST::one()/(2*ST::one());
1304 theta_ = theta_default_;
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);
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();
1330 myA(1,0) = as<Scalar>( one - theta_ );
1332 myb(0) = as<Scalar>( one - theta_ );
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);
1346 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList)
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_);
1353 this->setMyParamList(paramList);
1356 Scalar theta_default_;
1361 template<
class Scalar>
1362 class Implicit1Stage2ndOrderGauss_RKBT :
1363 virtual public RKButcherTableauDefaultBase<Scalar>
1366 Implicit1Stage2ndOrderGauss_RKBT()
1369 std::ostringstream myDescription;
1370 myDescription << Implicit1Stage2ndOrderGauss_name() <<
"\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" 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();
1395 this->setMyDescription(myDescription.str());
1399 this->setMy_order(2);
1404 template<
class Scalar>
1405 class Implicit2Stage4thOrderGauss_RKBT :
1406 virtual public RKButcherTableauDefaultBase<Scalar>
1409 Implicit2Stage4thOrderGauss_RKBT()
1412 std::ostringstream myDescription;
1413 myDescription << Implicit2Stage4thOrderGauss_name() <<
"\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;
1436 myA(0,0) = onefourth;
1437 myA(0,1) = onefourth-alpha;
1438 myA(1,0) = onefourth+alpha;
1439 myA(1,1) = onefourth;
1442 myc(0) = onehalf-alpha;
1443 myc(1) = onehalf+alpha;
1444 this->setMyDescription(myDescription.str());
1448 this->setMy_order(4);
1453 template<
class Scalar>
1454 class Implicit3Stage6thOrderGauss_RKBT :
1455 virtual public RKButcherTableauDefaultBase<Scalar>
1458 Implicit3Stage6thOrderGauss_RKBT()
1461 std::ostringstream myDescription;
1462 myDescription << Implicit3Stage6thOrderGauss_name() <<
"\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);
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());
1508 this->setMy_order(6);
1513 template<
class Scalar>
1514 class Implicit1Stage1stOrderRadauA_RKBT :
1515 virtual public RKButcherTableauDefaultBase<Scalar>
1518 Implicit1Stage1stOrderRadauA_RKBT()
1521 std::ostringstream myDescription;
1522 myDescription << Implicit1Stage1stOrderRadauA_name() <<
"\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" 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();
1542 this->setMyDescription(myDescription.str());
1546 this->setMy_order(1);
1551 template<
class Scalar>
1552 class Implicit2Stage3rdOrderRadauA_RKBT :
1553 virtual public RKButcherTableauDefaultBase<Scalar>
1556 Implicit2Stage3rdOrderRadauA_RKBT()
1559 std::ostringstream myDescription;
1560 myDescription << Implicit2Stage3rdOrderRadauA_name() <<
"\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));
1585 myc(1) = as<Scalar>(2*one/(3*one));
1586 this->setMyDescription(myDescription.str());
1590 this->setMy_order(3);
1595 template<
class Scalar>
1596 class Implicit3Stage5thOrderRadauA_RKBT :
1597 virtual public RKButcherTableauDefaultBase<Scalar>
1600 Implicit3Stage5thOrderRadauA_RKBT()
1603 std::ostringstream myDescription;
1604 myDescription << Implicit3Stage5thOrderRadauA_name() <<
"\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) );
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());
1642 this->setMy_order(5);
1647 template<
class Scalar>
1648 class Implicit1Stage1stOrderRadauB_RKBT :
1649 virtual public RKButcherTableauDefaultBase<Scalar>
1652 Implicit1Stage1stOrderRadauB_RKBT()
1655 std::ostringstream myDescription;
1656 myDescription << Implicit1Stage1stOrderRadauB_name() <<
"\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" 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();
1675 this->setMyDescription(myDescription.str());
1679 this->setMy_order(1);
1684 template<
class Scalar>
1685 class Implicit2Stage3rdOrderRadauB_RKBT :
1686 virtual public RKButcherTableauDefaultBase<Scalar>
1689 Implicit2Stage3rdOrderRadauB_RKBT()
1692 std::ostringstream myDescription;
1693 myDescription << Implicit2Stage3rdOrderRadauB_name() <<
"\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" 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) );
1718 this->setMyDescription(myDescription.str());
1722 this->setMy_order(3);
1727 template<
class Scalar>
1728 class Implicit3Stage5thOrderRadauB_RKBT :
1729 virtual public RKButcherTableauDefaultBase<Scalar>
1732 Implicit3Stage5thOrderRadauB_RKBT()
1735 std::ostringstream myDescription;
1736 myDescription << Implicit3Stage5thOrderRadauB_name() <<
"\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" 1754 <<
"b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'" 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) );
1777 this->setMyDescription(myDescription.str());
1781 this->setMy_order(5);
1786 template<
class Scalar>
1787 class Implicit2Stage2ndOrderLobattoA_RKBT :
1788 virtual public RKButcherTableauDefaultBase<Scalar>
1791 Implicit2Stage2ndOrderLobattoA_RKBT()
1794 std::ostringstream myDescription;
1795 myDescription << Implicit2Stage2ndOrderLobattoA_name() <<
"\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" 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();
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) );
1821 this->setMyDescription(myDescription.str());
1825 this->setMy_order(2);
1830 template<
class Scalar>
1831 class Implicit3Stage4thOrderLobattoA_RKBT :
1832 virtual public RKButcherTableauDefaultBase<Scalar>
1835 Implicit3Stage4thOrderLobattoA_RKBT()
1838 std::ostringstream myDescription;
1839 myDescription << Implicit3Stage4thOrderLobattoA_name() <<
"\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();
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) );
1871 myc(1) = as<Scalar>( one/(2*one) );
1873 this->setMyDescription(myDescription.str());
1877 this->setMy_order(4);
1882 template<
class Scalar>
1883 class Implicit4Stage6thOrderLobattoA_RKBT :
1884 virtual public RKButcherTableauDefaultBase<Scalar>
1887 Implicit4Stage6thOrderLobattoA_RKBT()
1891 typedef Teuchos::ScalarTraits<Scalar> ST;
1893 std::ostringstream myDescription;
1894 myDescription << Implicit4Stage6thOrderLobattoA_name() <<
"\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" 1904 <<
" [ (11+sqrt(5)/120 ]\n" 1905 <<
" [ (11-sqrt(5)/120 ]\n" 1908 <<
" [ (25-sqrt(5))/120 ]\n" 1909 <<
" [ (25+13*sqrt(5))/120 ]\n" 1912 <<
" [ (25-13*sqrt(5))/120 ]\n" 1913 <<
" [ (25+sqrt(5))/120 ]\n" 1916 <<
" [ (-1+sqrt(5))/120 ]\n" 1917 <<
" [ (-1-sqrt(5))/120 ]\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();
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) );
1948 myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
1949 myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
1951 this->setMyDescription(myDescription.str());
1955 this->setMy_order(6);
1960 template<
class Scalar>
1961 class Implicit2Stage2ndOrderLobattoB_RKBT :
1962 virtual public RKButcherTableauDefaultBase<Scalar>
1965 Implicit2Stage2ndOrderLobattoB_RKBT()
1968 std::ostringstream myDescription;
1969 myDescription << Implicit2Stage2ndOrderLobattoB_name() <<
"\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" 1977 <<
"A = [ 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) );
1989 myA(1,0) = as<Scalar>( one/(2*one) );
1991 myb(0) = as<Scalar>( one/(2*one) );
1992 myb(1) = as<Scalar>( one/(2*one) );
1995 this->setMyDescription(myDescription.str());
1999 this->setMy_order(2);
2004 template<
class Scalar>
2005 class Implicit3Stage4thOrderLobattoB_RKBT :
2006 virtual public RKButcherTableauDefaultBase<Scalar>
2009 Implicit3Stage4thOrderLobattoB_RKBT()
2012 std::ostringstream myDescription;
2013 myDescription << Implicit3Stage4thOrderLobattoB_name() <<
"\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) );
2035 myA(1,0) = as<Scalar>( one/(6*one) );
2036 myA(1,1) = as<Scalar>( one/(3*one) );
2038 myA(2,0) = as<Scalar>( one/(6*one) );
2039 myA(2,1) = as<Scalar>( 5*one/(6*one) );
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) );
2045 myc(1) = as<Scalar>( one/(2*one) );
2047 this->setMyDescription(myDescription.str());
2051 this->setMy_order(4);
2056 template<
class Scalar>
2057 class Implicit4Stage6thOrderLobattoB_RKBT :
2058 virtual public RKButcherTableauDefaultBase<Scalar>
2061 Implicit4Stage6thOrderLobattoB_RKBT()
2064 std::ostringstream myDescription;
2065 myDescription << Implicit4Stage6thOrderLobattoB_name() <<
"\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) );
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) );
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) );
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) );
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) );
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) );
2109 this->setMyDescription(myDescription.str());
2113 this->setMy_order(6);
2118 template<
class Scalar>
2119 class Implicit2Stage2ndOrderLobattoC_RKBT :
2120 virtual public RKButcherTableauDefaultBase<Scalar>
2123 Implicit2Stage2ndOrderLobattoC_RKBT()
2126 std::ostringstream myDescription;
2127 myDescription << Implicit2Stage2ndOrderLobattoC_name() <<
"\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" 2135 <<
"A = [ 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) );
2153 this->setMyDescription(myDescription.str());
2157 this->setMy_order(2);
2162 template<
class Scalar>
2163 class Implicit3Stage4thOrderLobattoC_RKBT :
2164 virtual public RKButcherTableauDefaultBase<Scalar>
2167 Implicit3Stage4thOrderLobattoC_RKBT()
2170 std::ostringstream myDescription;
2171 myDescription << Implicit3Stage4thOrderLobattoC_name() <<
"\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) );
2203 myc(1) = as<Scalar>( one/(2*one) );
2205 this->setMyDescription(myDescription.str());
2209 this->setMy_order(4);
2214 template<
class Scalar>
2215 class Implicit4Stage6thOrderLobattoC_RKBT :
2216 virtual public RKButcherTableauDefaultBase<Scalar>
2219 Implicit4Stage6thOrderLobattoC_RKBT()
2222 std::ostringstream myDescription;
2223 myDescription << Implicit4Stage6thOrderLobattoC_name() <<
"\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) );
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) );
2267 this->setMyDescription(myDescription.str());
2271 this->setMy_order(6);
2277 template<
class Scalar>
2278 class SDIRK5Stage5thOrder_RKBT :
2279 virtual public RKButcherTableauDefaultBase<Scalar>
2282 SDIRK5Stage5thOrder_RKBT()
2285 std::ostringstream myDescription;
2286 myDescription << SDIRK5Stage5thOrder_name() <<
"\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" 2293 <<
"c = [ (6-sqrt(6))/10 ]\n" 2294 <<
" [ (6+9*sqrt(6))/35 ]\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" 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" 2311 <<
" [ (6-sqrt(6))/10 ]\n" 2312 <<
" [ (-267+88*sqrt(6))/500 ]\n" 2313 <<
" [ (1329-544*sqrt(6))/2500 ]\n" 2317 <<
" [ (6-sqrt(6))/10 ]\n" 2318 <<
" [ (-96+131*sqrt(6))/625 ]\n" 2323 <<
" [ (6-sqrt(6))/10 ]\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) );
2344 myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
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) );
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) );
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) );
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) );
2375 myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
2377 myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
2378 myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
2380 this->setMyDescription(myDescription.str());
2384 this->setMy_order(5);
2389 template<
class Scalar>
2390 class SDIRK5Stage4thOrder_RKBT :
2391 virtual public RKButcherTableauDefaultBase<Scalar>
2394 SDIRK5Stage4thOrder_RKBT()
2397 std::ostringstream myDescription;
2398 myDescription << SDIRK5Stage4thOrder_name() <<
"\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" 2405 <<
"c = [ 1/4 3/4 11/20 1/2 1 ]'\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;
2427 myA(1,0) = as<Scalar>( one / (2*one) );
2428 myA(1,1) = onequarter;
2433 myA(2,0) = as<Scalar>( 17*one/(50*one) );
2434 myA(2,1) = as<Scalar>( -one/(25*one) );
2435 myA(2,2) = onequarter;
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;
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;
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;
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) );
2471 this->setMyDescription(myDescription.str());
2475 this->setMy_order(4);
2480 template<
class Scalar>
2481 class SDIRK3Stage4thOrder_RKBT :
2482 virtual public RKButcherTableauDefaultBase<Scalar>
2485 SDIRK3Stage4thOrder_RKBT()
2488 std::ostringstream myDescription;
2489 myDescription << SDIRK3Stage4thOrder_name() <<
"\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" 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)) );
2517 myA(1,0) = as<Scalar>( one/(2*one) - gamma );
2521 myA(2,0) = as<Scalar>( 2*gamma );
2522 myA(2,1) = as<Scalar>( one - 4*gamma );
2526 myb(1) = as<Scalar>( one-2*delta );
2530 myc(1) = as<Scalar>( one/(2*one) );
2531 myc(2) = as<Scalar>( one - gamma );
2533 this->setMyDescription(myDescription.str());
2537 this->setMy_order(4);
2545 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP