Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_RKButcherTableau.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 
00030 #ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HPP
00031 #define RYTHMOS_RK_BUTCHER_TABLEAU_HPP
00032 
00033 // disable clang warnings
00034 #ifdef __clang__
00035 #pragma clang system_header
00036 #endif
00037 
00038 #include "Rythmos_Types.hpp"
00039 #include "Rythmos_RKButcherTableauBase.hpp"
00040 
00041 #include "Teuchos_Assert.hpp"
00042 #include "Teuchos_as.hpp"
00043 #include "Teuchos_StandardParameterEntryValidators.hpp"
00044 #include "Teuchos_Describable.hpp"
00045 #include "Teuchos_VerboseObject.hpp"
00046 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00047 #include "Teuchos_ParameterListAcceptor.hpp"
00048 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00049 
00050 #include "Thyra_ProductVectorBase.hpp"
00051 
00052 namespace Rythmos {
00053 
00054   using Teuchos::as;
00055 
00056   inline const std::string RKBT_ForwardEuler_name() { return  "Forward Euler"; } // done
00057   inline const std::string RKBT_BackwardEuler_name() { return  "Backward Euler"; } // done
00058   inline const std::string Explicit4Stage_name() { return  "Explicit 4 Stage"; } // done
00059   inline const std::string Explicit3_8Rule_name() { return  "Explicit 3/8 Rule"; } // done
00060 
00061   inline const std::string ExplicitTrapezoidal_name() { return  "Explicit Trapezoidal"; } // done
00062   inline const std::string Explicit2Stage2ndOrderRunge_name() { return  "Explicit 2 Stage 2nd order by Runge"; } // done
00063   inline const std::string Explicit3Stage3rdOrderHeun_name() { return  "Explicit 3 Stage 3rd order by Heun"; } // done
00064   inline const std::string Explicit3Stage3rdOrder_name() { return  "Explicit 3 Stage 3rd order"; } // done
00065   inline const std::string Explicit3Stage3rdOrderTVD_name() { return  "Explicit 3 Stage 3rd order TVD"; } // done
00066   inline const std::string Explicit4Stage3rdOrderRunge_name() { return  "Explicit 4 Stage 3rd order by Runge"; } // done
00067 
00068   inline const std::string IRK1StageTheta_name() { return  "IRK 1 Stage Theta Method"; } // done
00069   inline const std::string IRK2StageTheta_name() { return  "IRK 2 Stage Theta Method"; } // done
00070   inline const std::string Implicit1Stage2ndOrderGauss_name() { return  "Implicit 1 Stage 2nd order Gauss"; } // done
00071   inline const std::string Implicit2Stage4thOrderGauss_name() { return  "Implicit 2 Stage 4th order Gauss"; } // done
00072   inline const std::string Implicit3Stage6thOrderGauss_name() { return  "Implicit 3 Stage 6th order Gauss"; } // done
00073 
00074   inline const std::string Implicit1Stage1stOrderRadauA_name() { return  "Implicit 1 Stage 1st order Radau left"; } // done
00075   inline const std::string Implicit2Stage3rdOrderRadauA_name() { return  "Implicit 2 Stage 3rd order Radau left"; } // done
00076   inline const std::string Implicit3Stage5thOrderRadauA_name() { return  "Implicit 3 Stage 5th order Radau left"; } // done
00077 
00078   inline const std::string Implicit1Stage1stOrderRadauB_name() { return  "Implicit 1 Stage 1st order Radau right"; } // done
00079   inline const std::string Implicit2Stage3rdOrderRadauB_name() { return  "Implicit 2 Stage 3rd order Radau right"; } // done
00080   inline const std::string Implicit3Stage5thOrderRadauB_name() { return  "Implicit 3 Stage 5th order Radau right"; } // done
00081 
00082   inline const std::string Implicit2Stage2ndOrderLobattoA_name() { return  "Implicit 2 Stage 2nd order Lobatto A"; } // done
00083   inline const std::string Implicit3Stage4thOrderLobattoA_name() { return  "Implicit 3 Stage 4th order Lobatto A"; } // done
00084   inline const std::string Implicit4Stage6thOrderLobattoA_name() { return  "Implicit 4 Stage 6th order Lobatto A"; } // done
00085 
00086   inline const std::string Implicit2Stage2ndOrderLobattoB_name() { return  "Implicit 2 Stage 2nd order Lobatto B"; } // done
00087   inline const std::string Implicit3Stage4thOrderLobattoB_name() { return  "Implicit 3 Stage 4th order Lobatto B"; } // done
00088   inline const std::string Implicit4Stage6thOrderLobattoB_name() { return  "Implicit 4 Stage 6th order Lobatto B"; } // done
00089 
00090   inline const std::string Implicit2Stage2ndOrderLobattoC_name() { return  "Implicit 2 Stage 2nd order Lobatto C"; } // done
00091   inline const std::string Implicit3Stage4thOrderLobattoC_name() { return  "Implicit 3 Stage 4th order Lobatto C"; } // done
00092   inline const std::string Implicit4Stage6thOrderLobattoC_name() { return  "Implicit 4 Stage 6th order Lobatto C"; } // done
00093 
00094   inline const std::string Implicit2Stage4thOrderHammerHollingsworth_name() { return  "Implicit 2 Stage 4th Order Hammer & Hollingsworth"; } // done
00095   inline const std::string Implicit3Stage6thOrderKuntzmannButcher_name() { return  "Implicit 3 Stage 6th Order Kuntzmann & Butcher"; } // done
00096   inline const std::string Implicit4Stage8thOrderKuntzmannButcher_name() { return  "Implicit 4 Stage 8th Order Kuntzmann & Butcher"; } // done
00097 
00098   inline const std::string DIRK2Stage3rdOrder_name() { return  "Diagonal IRK 2 Stage 3rd order"; } // done
00099 
00100   inline const std::string SDIRK2Stage2ndOrder_name() { return  "Singly Diagonal IRK 2 Stage 2nd order"; } // done
00101   inline const std::string SDIRK2Stage3rdOrder_name() { return  "Singly Diagonal IRK 2 Stage 3rd order"; } // done
00102   inline const std::string SDIRK5Stage5thOrder_name() { return  "Singly Diagonal IRK 5 Stage 5th order"; } // done
00103   inline const std::string SDIRK5Stage4thOrder_name() { return  "Singly Diagonal IRK 5 Stage 4th order"; } // done
00104   inline const std::string SDIRK3Stage4thOrder_name() { return  "Singly Diagonal IRK 3 Stage 4th order"; } // done
00105 
00106 template<class Scalar>
00107 class RKButcherTableauDefaultBase :
00108   virtual public RKButcherTableauBase<Scalar>,
00109   virtual public Teuchos::ParameterListAcceptorDefaultBase
00110 {
00111   public:
00113     virtual int numStages() const { return A_.numRows(); }
00115     virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const { return A_; }
00117     virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const { return b_; }
00119     virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const { return c_; }
00121     virtual int order() const { return order_; }
00123     virtual void setDescription(std::string longDescription) { longDescription_ = longDescription; }
00124 
00126     virtual void initialize(
00127         const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
00128         const Teuchos::SerialDenseVector<int,Scalar>& b_in,
00129         const Teuchos::SerialDenseVector<int,Scalar>& c_in,
00130         const int order_in,
00131         const std::string& longDescription_in
00132         )
00133     {
00134       const int numStages_in = A_in.numRows();
00135       TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
00136       TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
00137       TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
00138       TEUCHOS_ASSERT_EQUALITY( c_in.length(), numStages_in );
00139       TEUCHOS_ASSERT( order_in > 0 );
00140       A_ = A_in;
00141       b_ = b_in;
00142       c_ = c_in;
00143       order_ = order_in;
00144       longDescription_ = longDescription_in;
00145     }
00146 
00147     /* \brief Redefined from Teuchos::ParameterListAcceptorDefaultBase */
00149 
00151     virtual void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00152     {
00153       TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
00154       paramList->validateParameters(*this->getValidParameters());
00155       Teuchos::readVerboseObjectSublist(&*paramList,this);
00156       setMyParamList(paramList);
00157     }
00158 
00160     virtual RCP<const Teuchos::ParameterList> getValidParameters() const
00161     {
00162       if (is_null(validPL_)) {
00163         validPL_ = Teuchos::parameterList();
00164         validPL_->set("Description","",this->getMyDescription());
00165         Teuchos::setupVerboseObjectSublist(&*validPL_);
00166       }
00167       return validPL_;
00168     }
00169 
00171 
00172     /* \brief Redefined from Teuchos::Describable */
00174 
00176     virtual std::string description() const { return "Rythmos::RKButcherTableauDefaultBase"; }
00177 
00179     virtual void describe(
00180       Teuchos::FancyOStream &out,
00181       const Teuchos::EVerbosityLevel verbLevel
00182       ) const
00183     {
00184       if (verbLevel != Teuchos::VERB_NONE) {
00185         out << this->description() << std::endl;
00186         out << this->getMyDescription() << std::endl;
00187         out << "number of Stages = " << this->numStages() << std::endl;
00188         out << "A = " << this->A() << std::endl;
00189         out << "b = " << this->b() << std::endl;
00190         out << "c = " << this->c() << std::endl;
00191         out << "order = " << this->order() << std::endl;
00192       }
00193     }
00194 
00196 
00197   protected:
00198     void setMyDescription(std::string longDescription) { longDescription_ = longDescription; }
00199     const std::string& getMyDescription() const { return longDescription_; }
00200 
00201     void setMy_A(const Teuchos::SerialDenseMatrix<int,Scalar>& new_A) { A_ = new_A; }
00202     void setMy_b(const Teuchos::SerialDenseVector<int,Scalar>& new_b) { b_ = new_b; }
00203     void setMy_c(const Teuchos::SerialDenseVector<int,Scalar>& new_c) { c_ = new_c; }
00204     void setMy_order(const int& new_order) { order_ = new_order; }
00205 
00206     void setMyValidParameterList( const RCP<ParameterList> validPL ) { validPL_ = validPL; }
00207     RCP<ParameterList> getMyNonconstValidParameterList() { return validPL_; }
00208 
00209   private:
00210     Teuchos::SerialDenseMatrix<int,Scalar> A_;
00211     Teuchos::SerialDenseVector<int,Scalar> b_;
00212     Teuchos::SerialDenseVector<int,Scalar> c_;
00213     int order_;
00214     std::string longDescription_;
00215     mutable RCP<ParameterList> validPL_;
00216 };
00217 
00218 
00219 // Nonmember constructor
00220 template<class Scalar>
00221 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
00222 {
00223   return(rcp(new RKButcherTableauDefaultBase<Scalar>()));
00224 }
00225 
00226 // Nonmember constructor
00227 template<class Scalar>
00228 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau(
00229     const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
00230     const Teuchos::SerialDenseVector<int,Scalar>& b_in,
00231     const Teuchos::SerialDenseVector<int,Scalar>& c_in,
00232     int order_in,
00233     const std::string& description_in = ""
00234     )
00235 {
00236   RCP<RKButcherTableauDefaultBase<Scalar> > rkbt = rcp(new RKButcherTableauDefaultBase<Scalar>());
00237   rkbt->initialize(A_in,b_in,c_in,order_in,description_in);
00238   return(rkbt);
00239 }
00240 
00241 
00242 template<class Scalar>
00243 class BackwardEuler_RKBT :
00244   virtual public RKButcherTableauDefaultBase<Scalar>
00245 {
00246   public:
00247   BackwardEuler_RKBT()
00248   {
00249     std::ostringstream myDescription;
00250     myDescription << RKBT_BackwardEuler_name() << "\n"
00251                 << "c = [ 1 ]'\n"
00252                 << "A = [ 1 ]\n"
00253                 << "b = [ 1 ]'" << std::endl;
00254     typedef ScalarTraits<Scalar> ST;
00255     Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
00256     myA(0,0) = ST::one();
00257     Teuchos::SerialDenseVector<int,Scalar> myb(1);
00258     myb(0) = ST::one();
00259     Teuchos::SerialDenseVector<int,Scalar> myc(1);
00260     myc(0) = ST::one();
00261 
00262     this->setMyDescription(myDescription.str());
00263     this->setMy_A(myA);
00264     this->setMy_b(myb);
00265     this->setMy_c(myc);
00266     this->setMy_order(1);
00267   }
00268 };
00269 
00270 
00271 
00272 template<class Scalar>
00273 class ForwardEuler_RKBT :
00274   virtual public RKButcherTableauDefaultBase<Scalar>
00275 {
00276   public:
00277 
00278     ForwardEuler_RKBT()
00279     {
00280       std::ostringstream myDescription;
00281       myDescription << RKBT_ForwardEuler_name() << "\n"
00282                   << "c = [ 0 ]'\n"
00283                   << "A = [ 0 ]\n"
00284                   << "b = [ 1 ]'" << std::endl;
00285       typedef ScalarTraits<Scalar> ST;
00286       Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
00287       Teuchos::SerialDenseVector<int,Scalar> myb(1);
00288       myb(0) = ST::one();
00289       Teuchos::SerialDenseVector<int,Scalar> myc(1);
00290 
00291       this->setMyDescription(myDescription.str());
00292       this->setMy_A(myA);
00293       this->setMy_b(myb);
00294       this->setMy_c(myc);
00295       this->setMy_order(1);
00296     }
00297 };
00298 
00299 
00300 template<class Scalar>
00301 class Explicit4Stage4thOrder_RKBT :
00302   virtual public RKButcherTableauDefaultBase<Scalar>
00303 {
00304   public:
00305     Explicit4Stage4thOrder_RKBT()
00306     {
00307       std::ostringstream myDescription;
00308       myDescription << Explicit4Stage_name() << "\n"
00309                   << "\"The\" Runge-Kutta Method (explicit):\n"
00310                   << "Solving Ordinary Differential Equations I:\n"
00311                   << "Nonstiff Problems, 2nd Revised Edition\n"
00312                   << "E. Hairer, S.P. Norsett, G. Wanner\n"
00313                   << "Table 1.2, pg 138\n"
00314                   << "c = [  0  1/2 1/2  1  ]'\n"
00315                   << "A = [  0              ] \n"
00316                   << "    [ 1/2  0          ]\n"
00317                   << "    [  0  1/2  0      ]\n"
00318                   << "    [  0   0   1   0  ]\n"
00319                   << "b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
00320       typedef ScalarTraits<Scalar> ST;
00321       Scalar one = ST::one();
00322       Scalar zero = ST::zero();
00323       Scalar onehalf = ST::one()/(2*ST::one());
00324       Scalar onesixth = ST::one()/(6*ST::one());
00325       Scalar onethird = ST::one()/(3*ST::one());
00326 
00327       int myNumStages = 4;
00328       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00329       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00330       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00331 
00332       // Fill A:
00333       myA(0,0) = zero;
00334       myA(0,1) = zero;
00335       myA(0,2) = zero;
00336       myA(0,3) = zero;
00337 
00338       myA(1,0) = onehalf;
00339       myA(1,1) = zero;
00340       myA(1,2) = zero;
00341       myA(1,3) = zero;
00342 
00343       myA(2,0) = zero;
00344       myA(2,1) = onehalf;
00345       myA(2,2) = zero;
00346       myA(2,3) = zero;
00347 
00348       myA(3,0) = zero;
00349       myA(3,1) = zero;
00350       myA(3,2) = one;
00351       myA(3,3) = zero;
00352 
00353       // Fill myb:
00354       myb(0) = onesixth;
00355       myb(1) = onethird;
00356       myb(2) = onethird;
00357       myb(3) = onesixth;
00358 
00359       // fill b_c_
00360       myc(0) = zero;
00361       myc(1) = onehalf;
00362       myc(2) = onehalf;
00363       myc(3) = one;
00364 
00365       this->setMyDescription(myDescription.str());
00366       this->setMy_A(myA);
00367       this->setMy_b(myb);
00368       this->setMy_c(myc);
00369       this->setMy_order(4);
00370     }
00371 };
00372 
00373 
00374 template<class Scalar>
00375 class Explicit3_8Rule_RKBT :
00376   virtual public RKButcherTableauDefaultBase<Scalar>
00377 {
00378   public:
00379     Explicit3_8Rule_RKBT()
00380     {
00381 
00382       std::ostringstream myDescription;
00383       myDescription << Explicit3_8Rule_name() << "\n"
00384                   << "Solving Ordinary Differential Equations I:\n"
00385                   << "Nonstiff Problems, 2nd Revised Edition\n"
00386                   << "E. Hairer, S.P. Norsett, G. Wanner\n"
00387                   << "Table 1.2, pg 138\n"
00388                   << "c = [  0  1/3 2/3  1  ]'\n"
00389                   << "A = [  0              ]\n"
00390                   << "    [ 1/3  0          ]\n"
00391                   << "    [-1/3  1   0      ]\n"
00392                   << "    [  1  -1   1   0  ]\n"
00393                   << "b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
00394       typedef ScalarTraits<Scalar> ST;
00395       int myNumStages = 4;
00396       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00397       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00398       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00399 
00400       Scalar one = ST::one();
00401       Scalar zero = ST::zero();
00402       Scalar one_third    = as<Scalar>(ST::one()/(3*ST::one()));
00403       Scalar two_third    = as<Scalar>(2*ST::one()/(3*ST::one()));
00404       Scalar one_eighth   = as<Scalar>(ST::one()/(8*ST::one()));
00405       Scalar three_eighth = as<Scalar>(3*ST::one()/(8*ST::one()));
00406 
00407       // Fill myA:
00408       myA(0,0) = zero;
00409       myA(0,1) = zero;
00410       myA(0,2) = zero;
00411       myA(0,3) = zero;
00412 
00413       myA(1,0) = one_third;
00414       myA(1,1) = zero;
00415       myA(1,2) = zero;
00416       myA(1,3) = zero;
00417 
00418       myA(2,0) = as<Scalar>(-one_third);
00419       myA(2,1) = one;
00420       myA(2,2) = zero;
00421       myA(2,3) = zero;
00422 
00423       myA(3,0) = one;
00424       myA(3,1) = as<Scalar>(-one);
00425       myA(3,2) = one;
00426       myA(3,3) = zero;
00427 
00428       // Fill myb:
00429       myb(0) = one_eighth;
00430       myb(1) = three_eighth;
00431       myb(2) = three_eighth;
00432       myb(3) = one_eighth;
00433 
00434       // Fill myc:
00435       myc(0) = zero;
00436       myc(1) = one_third;
00437       myc(2) = two_third;
00438       myc(3) = one;
00439 
00440       this->setMyDescription(myDescription.str());
00441       this->setMy_A(myA);
00442       this->setMy_b(myb);
00443       this->setMy_c(myc);
00444       this->setMy_order(4);
00445     }
00446 };
00447 
00448 
00449 template<class Scalar>
00450 class Explicit4Stage3rdOrderRunge_RKBT :
00451   virtual public RKButcherTableauDefaultBase<Scalar>
00452 {
00453   public:
00454     Explicit4Stage3rdOrderRunge_RKBT()
00455     {
00456 
00457       std::ostringstream myDescription;
00458       myDescription << Explicit4Stage3rdOrderRunge_name() << "\n"
00459                   << "Solving Ordinary Differential Equations I:\n"
00460                   << "Nonstiff Problems, 2nd Revised Edition\n"
00461                   << "E. Hairer, S.P. Norsett, G. Wanner\n"
00462                   << "Table 1.1, pg 135\n"
00463                   << "c = [  0  1/2  1   1  ]'\n"
00464                   << "A = [  0              ]\n"
00465                   << "    [ 1/2  0          ]\n"
00466                   << "    [  0   1   0      ]\n"
00467                   << "    [  0   0   1   0  ]\n"
00468                   << "b = [ 1/6 2/3  0  1/6 ]'" << std::endl;
00469       typedef ScalarTraits<Scalar> ST;
00470       int myNumStages = 4;
00471       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00472       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00473       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00474 
00475       Scalar one = ST::one();
00476       Scalar onehalf = ST::one()/(2*ST::one());
00477       Scalar onesixth = ST::one()/(6*ST::one());
00478       Scalar twothirds = 2*ST::one()/(3*ST::one());
00479       Scalar zero = ST::zero();
00480 
00481       // Fill A:
00482       myA(0,0) = zero;
00483       myA(0,1) = zero;
00484       myA(0,2) = zero;
00485       myA(0,3) = zero;
00486 
00487       myA(1,0) = onehalf;
00488       myA(1,1) = zero;
00489       myA(1,2) = zero;
00490       myA(1,3) = zero;
00491 
00492       myA(2,0) = zero;
00493       myA(2,1) = one;
00494       myA(2,2) = zero;
00495       myA(2,3) = zero;
00496 
00497       myA(3,0) = zero;
00498       myA(3,1) = zero;
00499       myA(3,2) = one;
00500       myA(3,3) = zero;
00501 
00502       // Fill b:
00503       myb(0) = onesixth;
00504       myb(1) = twothirds;
00505       myb(2) = zero;
00506       myb(3) = onesixth;
00507 
00508       // Fill myc:
00509       myc(0) = zero;
00510       myc(1) = onehalf;
00511       myc(2) = one;
00512       myc(3) = one;
00513 
00514       this->setMyDescription(myDescription.str());
00515       this->setMy_A(myA);
00516       this->setMy_b(myb);
00517       this->setMy_c(myc);
00518       this->setMy_order(3);
00519     }
00520 };
00521 
00522 
00523 template<class Scalar>
00524 class Explicit3Stage3rdOrder_RKBT :
00525   virtual public RKButcherTableauDefaultBase<Scalar>
00526 {
00527   public:
00528     Explicit3Stage3rdOrder_RKBT()
00529     {
00530 
00531       std::ostringstream myDescription;
00532       myDescription << Explicit3Stage3rdOrder_name() << "\n"
00533                   << "c = [  0  1/2  1  ]'\n"
00534                   << "A = [  0          ]\n"
00535                   << "    [ 1/2  0      ]\n"
00536                   << "    [ -1   2   0  ]\n"
00537                   << "b = [ 1/6 4/6 1/6 ]'" << std::endl;
00538       typedef ScalarTraits<Scalar> ST;
00539       Scalar one = ST::one();
00540       Scalar two = as<Scalar>(2*ST::one());
00541       Scalar zero = ST::zero();
00542       Scalar onehalf = ST::one()/(2*ST::one());
00543       Scalar onesixth = ST::one()/(6*ST::one());
00544       Scalar foursixth = 4*ST::one()/(6*ST::one());
00545 
00546       int myNumStages = 3;
00547       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00548       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00549       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00550 
00551       // Fill myA:
00552       myA(0,0) = zero;
00553       myA(0,1) = zero;
00554       myA(0,2) = zero;
00555 
00556       myA(1,0) = onehalf;
00557       myA(1,1) = zero;
00558       myA(1,2) = zero;
00559 
00560       myA(2,0) = -one;
00561       myA(2,1) = two;
00562       myA(2,2) = zero;
00563 
00564       // Fill myb:
00565       myb(0) = onesixth;
00566       myb(1) = foursixth;
00567       myb(2) = onesixth;
00568 
00569       // fill b_c_
00570       myc(0) = zero;
00571       myc(1) = onehalf;
00572       myc(2) = one;
00573 
00574       this->setMyDescription(myDescription.str());
00575       this->setMy_A(myA);
00576       this->setMy_b(myb);
00577       this->setMy_c(myc);
00578       this->setMy_order(3);
00579     }
00580 };
00581 
00582 
00583 template<class Scalar>
00584 class Explicit3Stage3rdOrderTVD_RKBT :
00585   virtual public RKButcherTableauDefaultBase<Scalar>
00586 {
00587   public:
00588     Explicit3Stage3rdOrderTVD_RKBT()
00589     {
00590 
00591       std::ostringstream myDescription;
00592       myDescription << Explicit3Stage3rdOrderTVD_name() << "\n"
00593                     << "Sigal Gottlieb and Chi-Wang Shu\n"
00594                     << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
00595                     << "Mathematics of Computation\n"
00596                     << "Volume 67, Number 221, January 1998, pp. 73-85\n"
00597                     << "c = [  0   1  1/2 ]'\n"
00598                     << "A = [  0          ]\n"
00599                     << "    [  1   0      ]\n"
00600                     << "    [ 1/4 1/4  0  ]\n"
00601                     << "b = [ 1/6 1/6 4/6 ]'\n"
00602                     << "This is also written in the following set of updates.\n"
00603                     << "u1 = u^n + dt L(u^n)\n"
00604                     << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
00605                     << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3"
00606                     << std::endl;
00607       typedef ScalarTraits<Scalar> ST;
00608       Scalar one = ST::one();
00609       Scalar zero = ST::zero();
00610       Scalar onehalf = ST::one()/(2*ST::one());
00611       Scalar onefourth = ST::one()/(4*ST::one());
00612       Scalar onesixth = ST::one()/(6*ST::one());
00613       Scalar foursixth = 4*ST::one()/(6*ST::one());
00614 
00615       int myNumStages = 3;
00616       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00617       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00618       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00619 
00620       // Fill myA:
00621       myA(0,0) = zero;
00622       myA(0,1) = zero;
00623       myA(0,2) = zero;
00624 
00625       myA(1,0) = one;
00626       myA(1,1) = zero;
00627       myA(1,2) = zero;
00628 
00629       myA(2,0) = onefourth;
00630       myA(2,1) = onefourth;
00631       myA(2,2) = zero;
00632 
00633       // Fill myb:
00634       myb(0) = onesixth;
00635       myb(1) = onesixth;
00636       myb(2) = foursixth;
00637 
00638       // fill b_c_
00639       myc(0) = zero;
00640       myc(1) = one;
00641       myc(2) = onehalf;
00642 
00643       this->setMyDescription(myDescription.str());
00644       this->setMy_A(myA);
00645       this->setMy_b(myb);
00646       this->setMy_c(myc);
00647       this->setMy_order(3);
00648     }
00649 };
00650 
00651 
00652 template<class Scalar>
00653 class Explicit3Stage3rdOrderHeun_RKBT :
00654   virtual public RKButcherTableauDefaultBase<Scalar>
00655 {
00656   public:
00657     Explicit3Stage3rdOrderHeun_RKBT()
00658     {
00659       std::ostringstream myDescription;
00660       myDescription << Explicit3Stage3rdOrderHeun_name() << "\n"
00661                   << "Solving Ordinary Differential Equations I:\n"
00662                   << "Nonstiff Problems, 2nd Revised Edition\n"
00663                   << "E. Hairer, S.P. Norsett, G. Wanner\n"
00664                   << "Table 1.1, pg 135\n"
00665                   << "c = [  0  1/3 2/3 ]'\n"
00666                   << "A = [  0          ] \n"
00667                   << "    [ 1/3  0      ]\n"
00668                   << "    [  0  2/3  0  ]\n"
00669                   << "b = [ 1/4  0  3/4 ]'" << std::endl;
00670       typedef ScalarTraits<Scalar> ST;
00671       Scalar one = ST::one();
00672       Scalar zero = ST::zero();
00673       Scalar onethird = one/(3*one);
00674       Scalar twothirds = 2*one/(3*one);
00675       Scalar onefourth = one/(4*one);
00676       Scalar threefourths = 3*one/(4*one);
00677 
00678       int myNumStages = 3;
00679       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00680       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00681       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00682 
00683       // Fill myA:
00684       myA(0,0) = zero;
00685       myA(0,1) = zero;
00686       myA(0,2) = zero;
00687 
00688       myA(1,0) = onethird;
00689       myA(1,1) = zero;
00690       myA(1,2) = zero;
00691 
00692       myA(2,0) = zero;
00693       myA(2,1) = twothirds;
00694       myA(2,2) = zero;
00695 
00696       // Fill myb:
00697       myb(0) = onefourth;
00698       myb(1) = zero;
00699       myb(2) = threefourths;
00700 
00701       // fill b_c_
00702       myc(0) = zero;
00703       myc(1) = onethird;
00704       myc(2) = twothirds;
00705 
00706       this->setMyDescription(myDescription.str());
00707       this->setMy_A(myA);
00708       this->setMy_b(myb);
00709       this->setMy_c(myc);
00710       this->setMy_order(3);
00711     }
00712 };
00713 
00714 
00715 template<class Scalar>
00716 class Explicit2Stage2ndOrderRunge_RKBT :
00717   virtual public RKButcherTableauDefaultBase<Scalar>
00718 {
00719   public:
00720     Explicit2Stage2ndOrderRunge_RKBT()
00721     {
00722       std::ostringstream myDescription;
00723       myDescription << Explicit2Stage2ndOrderRunge_name() << "\n"
00724                   << "Also known as Explicit Midpoint\n"
00725                   << "Solving Ordinary Differential Equations I:\n"
00726                   << "Nonstiff Problems, 2nd Revised Edition\n"
00727                   << "E. Hairer, S.P. Norsett, G. Wanner\n"
00728                   << "Table 1.1, pg 135\n"
00729                   << "c = [  0  1/2 ]'\n"
00730                   << "A = [  0      ]\n"
00731                   << "    [ 1/2  0  ]\n"
00732                   << "b = [  0   1  ]'" << std::endl;
00733       typedef ScalarTraits<Scalar> ST;
00734       Scalar one = ST::one();
00735       Scalar zero = ST::zero();
00736       Scalar onehalf = ST::one()/(2*ST::one());
00737 
00738       int myNumStages = 2;
00739       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00740       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00741       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00742 
00743       // Fill myA:
00744       myA(0,0) = zero;
00745       myA(0,1) = zero;
00746 
00747       myA(1,0) = onehalf;
00748       myA(1,1) = zero;
00749 
00750       // Fill myb:
00751       myb(0) = zero;
00752       myb(1) = one;
00753 
00754       // fill b_c_
00755       myc(0) = zero;
00756       myc(1) = onehalf;
00757 
00758       this->setMyDescription(myDescription.str());
00759       this->setMy_A(myA);
00760       this->setMy_b(myb);
00761       this->setMy_c(myc);
00762       this->setMy_order(2);
00763     }
00764 };
00765 
00766 
00767 template<class Scalar>
00768 class ExplicitTrapezoidal_RKBT :
00769   virtual public RKButcherTableauDefaultBase<Scalar>
00770 {
00771   public:
00772     ExplicitTrapezoidal_RKBT()
00773     {
00774       std::ostringstream myDescription;
00775       myDescription << ExplicitTrapezoidal_name() << "\n"
00776                   << "c = [  0   1  ]'\n"
00777                   << "A = [  0      ]\n"
00778                   << "    [  1   0  ]\n"
00779                   << "b = [ 1/2 1/2 ]'" << std::endl;
00780       typedef ScalarTraits<Scalar> ST;
00781       Scalar one = ST::one();
00782       Scalar zero = ST::zero();
00783       Scalar onehalf = ST::one()/(2*ST::one());
00784 
00785       int myNumStages = 2;
00786       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00787       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00788       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00789 
00790       // Fill myA:
00791       myA(0,0) = zero;
00792       myA(0,1) = zero;
00793 
00794       myA(1,0) = one;
00795       myA(1,1) = zero;
00796 
00797       // Fill myb:
00798       myb(0) = onehalf;
00799       myb(1) = onehalf;
00800 
00801       // fill b_c_
00802       myc(0) = zero;
00803       myc(1) = one;
00804 
00805       this->setMyDescription(myDescription.str());
00806       this->setMy_A(myA);
00807       this->setMy_b(myb);
00808       this->setMy_c(myc);
00809       this->setMy_order(2);
00810     }
00811 };
00812 
00813 
00814 template<class Scalar>
00815 class SDIRK2Stage2ndOrder_RKBT :
00816   virtual public RKButcherTableauDefaultBase<Scalar>
00817 {
00818   public:
00819     SDIRK2Stage2ndOrder_RKBT()
00820     {
00821       std::ostringstream myDescription;
00822       myDescription << SDIRK2Stage2ndOrder_name() << "\n"
00823                   << "Computer Methods for ODEs and DAEs\n"
00824                   << "U. M. Ascher and L. R. Petzold\n"
00825                   << "p. 106\n"
00826                   << "gamma = (2+-sqrt(2))/2\n"
00827                   << "c = [  gamma   1     ]'\n"
00828                   << "A = [  gamma   0     ]\n"
00829                   << "    [ 1-gamma  gamma ]\n"
00830                   << "b = [ 1-gamma  gamma ]'" << std::endl;
00831 
00832       this->setMyDescription(myDescription.str());
00833       typedef ScalarTraits<Scalar> ST;
00834       Scalar one = ST::one();
00835       gamma_default_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
00836       gamma_ = gamma_default_;
00837       this->setupData();
00838 
00839       RCP<ParameterList> validPL = Teuchos::parameterList();
00840       validPL->set("Description","",this->getMyDescription());
00841       validPL->set<double>("gamma",gamma_default_,
00842         "The default value is gamma = (2-sqrt(2))/2. "
00843         "This will produce an L-stable 2nd order method with the stage "
00844         "times within the timestep.  Other values of gamma will still "
00845         "produce an L-stable scheme, but will only be 1st order accurate.");
00846       Teuchos::setupVerboseObjectSublist(&*validPL);
00847       this->setMyValidParameterList(validPL);
00848     }
00849     void setupData()
00850     {
00851       typedef ScalarTraits<Scalar> ST;
00852       int myNumStages = 2;
00853       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00854       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00855       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00856       Scalar one = ST::one();
00857       Scalar zero = ST::zero();
00858       myA(0,0) = gamma_;
00859       myA(0,1) = zero;
00860       myA(1,0) = as<Scalar>( one - gamma_ );
00861       myA(1,1) = gamma_;
00862       myb(0) = as<Scalar>( one - gamma_ );
00863       myb(1) = gamma_;
00864       myc(0) = gamma_;
00865       myc(1) = one;
00866 
00867       this->setMy_A(myA);
00868       this->setMy_b(myb);
00869       this->setMy_c(myc);
00870       this->setMy_order(2);
00871     }
00872     void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00873     {
00874       TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
00875       paramList->validateParameters(*this->getValidParameters());
00876       Teuchos::readVerboseObjectSublist(&*paramList,this);
00877       gamma_ = paramList->get<double>("gamma",gamma_default_);
00878       this->setupData();
00879       this->setMyParamList(paramList);
00880     }
00881   private:
00882     Scalar gamma_default_;
00883     Scalar gamma_;
00884 };
00885 
00886 
00887 // 04/07/09 tscoffe:  I verified manually that the Convergence Testing passes
00888 // with gamma_default_ = -1.
00889 template<class Scalar>
00890 class SDIRK2Stage3rdOrder_RKBT :
00891   virtual public RKButcherTableauDefaultBase<Scalar>
00892 {
00893   public:
00894     SDIRK2Stage3rdOrder_RKBT()
00895     {
00896       std::ostringstream myDescription;
00897       myDescription << SDIRK2Stage3rdOrder_name() << "\n"
00898                   << "Solving Ordinary Differential Equations I:\n"
00899                   << "Nonstiff Problems, 2nd Revised Edition\n"
00900                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
00901                   << "Table 7.2, pg 207\n"
00902                   << "gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n"
00903                   << "gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n"
00904                   << "c = [  gamma     1-gamma  ]'\n"
00905                   << "A = [  gamma     0        ]\n"
00906                   << "    [ 1-2*gamma  gamma    ]\n"
00907                   << "b = [ 1/2        1/2      ]'" << std::endl;
00908 
00909       this->setMyDescription(myDescription.str());
00910       thirdOrderAStable_default_ = true;
00911       secondOrderLStable_default_ = false;
00912       thirdOrderAStable_ = true;
00913       secondOrderLStable_ = false;
00914       typedef ScalarTraits<Scalar> ST;
00915       Scalar one = ST::one();
00916       gamma_default_ = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
00917       gamma_ = gamma_default_;
00918       this->setupData();
00919 
00920       RCP<ParameterList> validPL = Teuchos::parameterList();
00921       validPL->set("Description","",this->getMyDescription());
00922       validPL->set("3rd Order A-stable",thirdOrderAStable_default_,
00923         "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain "
00924         "a 3rd order A-stable scheme. '3rd Order A-stable' and "
00925         "'2nd Order L-stable' can not both be true.");
00926       validPL->set("2nd Order L-stable",secondOrderLStable_default_,
00927         "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain "
00928         "a 2nd order L-stable scheme. '3rd Order A-stable' and "
00929         "'2nd Order L-stable' can not both be true.");
00930       validPL->set("gamma",gamma_default_,
00931         "If both '3rd Order A-stable' and '2nd Order L-stable' "
00932         "are false, gamma will be used. The default value is the "
00933         "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
00934 
00935       Teuchos::setupVerboseObjectSublist(&*validPL);
00936       this->setMyValidParameterList(validPL);
00937     }
00938     void setupData()
00939     {
00940       typedef ScalarTraits<Scalar> ST;
00941       int myNumStages = 2;
00942       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
00943       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
00944       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
00945       Scalar one = ST::one();
00946       Scalar zero = ST::zero();
00947       Scalar gamma;
00948       if (thirdOrderAStable_)
00949         gamma = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
00950       else if (secondOrderLStable_)
00951         gamma = as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
00952       else
00953         gamma = gamma_;
00954       myA(0,0) = gamma;
00955       myA(0,1) = zero;
00956       myA(1,0) = as<Scalar>( one - 2*gamma );
00957       myA(1,1) = gamma;
00958       myb(0) = as<Scalar>( one/(2*one) );
00959       myb(1) = as<Scalar>( one/(2*one) );
00960       myc(0) = gamma;
00961       myc(1) = as<Scalar>( one - gamma );
00962 
00963       this->setMy_A(myA);
00964       this->setMy_b(myb);
00965       this->setMy_c(myc);
00966       this->setMy_order(3);
00967     }
00968     void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00969     {
00970       TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
00971       paramList->validateParameters(*this->getValidParameters());
00972       Teuchos::readVerboseObjectSublist(&*paramList,this);
00973       thirdOrderAStable_  = paramList->get("3rd Order A-stable",
00974                                            thirdOrderAStable_default_);
00975       secondOrderLStable_ = paramList->get("2nd Order L-stable",
00976                                            secondOrderLStable_default_);
00977       TEUCHOS_TEST_FOR_EXCEPTION(
00978         thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
00979         "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
00980       gamma_ = paramList->get("gamma",gamma_default_);
00981 
00982       this->setupData();
00983       this->setMyParamList(paramList);
00984     }
00985   private:
00986     bool thirdOrderAStable_default_;
00987     bool thirdOrderAStable_;
00988     bool secondOrderLStable_default_;
00989     bool secondOrderLStable_;
00990     Scalar gamma_default_;
00991     Scalar gamma_;
00992 };
00993 
00994 
00995 template<class Scalar>
00996 class DIRK2Stage3rdOrder_RKBT :
00997   virtual public RKButcherTableauDefaultBase<Scalar>
00998 {
00999   public:
01000     DIRK2Stage3rdOrder_RKBT()
01001     {
01002 
01003       std::ostringstream myDescription;
01004       myDescription << DIRK2Stage3rdOrder_name() << "\n"
01005                   << "Hammer & Hollingsworth method\n"
01006                   << "Solving Ordinary Differential Equations I:\n"
01007                   << "Nonstiff Problems, 2nd Revised Edition\n"
01008                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
01009                   << "Table 7.1, pg 205\n"
01010                   << "c = [  0   2/3 ]'\n"
01011                   << "A = [  0    0  ]\n"
01012                   << "    [ 1/3  1/3 ]\n"
01013                   << "b = [ 1/4  3/4 ]'" << std::endl;
01014       typedef ScalarTraits<Scalar> ST;
01015       int myNumStages = 2;
01016       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01017       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01018       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01019       Scalar one = ST::one();
01020       Scalar zero = ST::zero();
01021       myA(0,0) = zero;
01022       myA(0,1) = zero;
01023       myA(1,0) = as<Scalar>( one/(3*one) );
01024       myA(1,1) = as<Scalar>( one/(3*one) );
01025       myb(0) = as<Scalar>( one/(4*one) );
01026       myb(1) = as<Scalar>( 3*one/(4*one) );
01027       myc(0) = zero;
01028       myc(1) = as<Scalar>( 2*one/(3*one) );
01029       this->setMyDescription(myDescription.str());
01030       this->setMy_A(myA);
01031       this->setMy_b(myb);
01032       this->setMy_c(myc);
01033       this->setMy_order(3);
01034     }
01035 };
01036 
01037 
01038 template<class Scalar>
01039 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
01040   virtual public RKButcherTableauDefaultBase<Scalar>
01041 {
01042   public:
01043     Implicit3Stage6thOrderKuntzmannButcher_RKBT()
01044     {
01045 
01046       std::ostringstream myDescription;
01047       myDescription << Implicit3Stage6thOrderKuntzmannButcher_name() << "\n"
01048                   << "Kuntzmann & Butcher method\n"
01049                   << "Solving Ordinary Differential Equations I:\n"
01050                   << "Nonstiff Problems, 2nd Revised Edition\n"
01051                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
01052                   << "Table 7.4, pg 209\n"
01053                   << "c = [ 1/2-sqrt(15)/10   1/2              1/2-sqrt(15)/10  ]'\n"
01054                   << "A = [ 5/36              2/9-sqrt(15)/15  5/36-sqrt(15)/30 ]\n"
01055                   << "    [ 5/36+sqrt(15)/24  2/9              5/36-sqrt(15)/24 ]\n"
01056                   << "    [ 5/36+sqrt(15)/30  2/9+sqrt(15)/15  5/36             ]\n"
01057                   << "b = [ 5/18              4/9              5/18             ]'" << std::endl;
01058       typedef ScalarTraits<Scalar> ST;
01059       int myNumStages = 3;
01060       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01061       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01062       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01063       Scalar one = ST::one();
01064       myA(0,0) = as<Scalar>( 5*one/(36*one) );
01065       myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
01066       myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
01067       myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
01068       myA(1,1) = as<Scalar>( 2*one/(9*one) );
01069       myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
01070       myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
01071       myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
01072       myA(2,2) = as<Scalar>( 5*one/(36*one) );
01073       myb(0) = as<Scalar>( 5*one/(18*one) );
01074       myb(1) = as<Scalar>( 4*one/(9*one) );
01075       myb(2) = as<Scalar>( 5*one/(18*one) );
01076       myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
01077       myc(1) = as<Scalar>( one/(2*one) );
01078       myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
01079       this->setMyDescription(myDescription.str());
01080       this->setMy_A(myA);
01081       this->setMy_b(myb);
01082       this->setMy_c(myc);
01083       this->setMy_order(6);
01084     }
01085 };
01086 
01087 
01088 template<class Scalar>
01089 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
01090   virtual public RKButcherTableauDefaultBase<Scalar>
01091 {
01092   public:
01093     Implicit4Stage8thOrderKuntzmannButcher_RKBT()
01094     {
01095 
01096       std::ostringstream myDescription;
01097       myDescription << Implicit4Stage8thOrderKuntzmannButcher_name() << "\n"
01098                   << "Kuntzmann & Butcher method\n"
01099                   << "Solving Ordinary Differential Equations I:\n"
01100                   << "Nonstiff Problems, 2nd Revised Edition\n"
01101                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
01102                   << "Table 7.5, pg 209\n"
01103                   << "c = [ 1/2-w2     1/2-w2p     1/2+w2p     1/2+w2    ]'\n"
01104                   << "A = [ w1         w1p-w3+w4p  w1p-w3-w4p  w1-w5     ]\n"
01105                   << "    [ w1-w3p+w4  w1p         w1p-w5p     w1-w3p-w4 ]\n"
01106                   << "    [ w1+w3p+w4  w1p+w5p     w1p         w1+w3p-w4 ]\n"
01107                   << "    [ w1+w5      w1p+w3+w4p  w1p+w3-w4p  w1        ]\n"
01108                   << "b = [ 2*w1       2*w1p       2*w1p       2*w1      ]'\n"
01109                   << "w1 = 1/8-sqrt(30)/144\n"
01110                   << "w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n"
01111                   << "w3 = w2*(1/6+sqrt(30)/24)\n"
01112                   << "w4 = w2*(1/21+5*sqrt(30)/168)\n"
01113                   << "w5 = w2-2*w3\n"
01114                   << "w1p = 1/8+sqrt(30)/144\n"
01115                   << "w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n"
01116                   << "w3p = w2*(1/6-sqrt(30)/24)\n"
01117                   << "w4p = w2*(1/21-5*sqrt(30)/168)\n"
01118                   << "w5p = w2p-2*w3p" << std::endl;
01119       typedef ScalarTraits<Scalar> ST;
01120       int myNumStages = 4;
01121       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01122       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01123       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01124       Scalar one = ST::one();
01125       Scalar onehalf = as<Scalar>( one/(2*one) );
01126       Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
01127       Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
01128       Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
01129       Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
01130       Scalar w5 = as<Scalar>( w2-2*w3 );
01131       Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
01132       Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
01133       Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
01134       Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
01135       Scalar w5p = as<Scalar>( w2p-2*w3p );
01136       myA(0,0) = w1;
01137       myA(0,1) = w1p-w3+w4p;
01138       myA(0,2) = w1p-w3-w4p;
01139       myA(0,3) = w1-w5;
01140       myA(1,0) = w1-w3p+w4;
01141       myA(1,1) = w1p;
01142       myA(1,2) = w1p-w5p;
01143       myA(1,3) = w1-w3p-w4;
01144       myA(2,0) = w1+w3p+w4;
01145       myA(2,1) = w1p+w5p;
01146       myA(2,2) = w1p;
01147       myA(2,3) = w1+w3p-w4;
01148       myA(3,0) = w1+w5;
01149       myA(3,1) = w1p+w3+w4p;
01150       myA(3,2) = w1p+w3-w4p;
01151       myA(3,3) = w1;
01152       myb(0) = 2*w1;
01153       myb(1) = 2*w1p;
01154       myb(2) = 2*w1p;
01155       myb(3) = 2*w1;
01156       myc(0) = onehalf - w2;
01157       myc(1) = onehalf - w2p;
01158       myc(2) = onehalf + w2p;
01159       myc(3) = onehalf + w2;
01160       this->setMyDescription(myDescription.str());
01161       this->setMy_A(myA);
01162       this->setMy_b(myb);
01163       this->setMy_c(myc);
01164       this->setMy_order(8);
01165     }
01166 };
01167 
01168 
01169 template<class Scalar>
01170 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
01171   virtual public RKButcherTableauDefaultBase<Scalar>
01172 {
01173   public:
01174     Implicit2Stage4thOrderHammerHollingsworth_RKBT()
01175     {
01176 
01177       std::ostringstream myDescription;
01178       myDescription << Implicit2Stage4thOrderHammerHollingsworth_name() << "\n"
01179                   << "Hammer & Hollingsworth method\n"
01180                   << "Solving Ordinary Differential Equations I:\n"
01181                   << "Nonstiff Problems, 2nd Revised Edition\n"
01182                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
01183                   << "Table 7.3, pg 207\n"
01184                   << "c = [ 1/2-sqrt(3)/6  1/2+sqrt(3)/6 ]'\n"
01185                   << "A = [ 1/4            1/4-sqrt(3)/6 ]\n"
01186                   << "    [ 1/4+sqrt(3)/6  1/4           ]\n"
01187                   << "b = [ 1/2            1/2           ]'" << std::endl;
01188       typedef ScalarTraits<Scalar> ST;
01189       int myNumStages = 2;
01190       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01191       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01192       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01193       Scalar one = ST::one();
01194       Scalar onequarter = as<Scalar>( one/(4*one) );
01195       Scalar onehalf = as<Scalar>( one/(2*one) );
01196       myA(0,0) = onequarter;
01197       myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
01198       myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
01199       myA(1,1) = onequarter;
01200       myb(0) = onehalf;
01201       myb(1) = onehalf;
01202       myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
01203       myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
01204       this->setMyDescription(myDescription.str());
01205       this->setMy_A(myA);
01206       this->setMy_b(myb);
01207       this->setMy_c(myc);
01208       this->setMy_order(4);
01209     }
01210 };
01211 
01212 
01213 template<class Scalar>
01214 class IRK1StageTheta_RKBT :
01215   virtual public RKButcherTableauDefaultBase<Scalar>
01216 {
01217   public:
01218     IRK1StageTheta_RKBT()
01219     {
01220 
01221       std::ostringstream myDescription;
01222       myDescription << IRK1StageTheta_name() << "\n"
01223                   << "Non-standard finite-difference methods\n"
01224                   << "in dynamical systems, P. Kama,\n"
01225                   << "Dissertation, University of Pretoria, pg. 49.\n"
01226                   << "Comment:  Generalized Implicit Midpoint Method\n"
01227                   << "c = [ theta ]'\n"
01228                   << "A = [ theta ]\n"
01229                   << "b = [  1  ]'\n"
01230                   << std::endl;
01231 
01232       this->setMyDescription(myDescription.str());
01233       typedef ScalarTraits<Scalar> ST;
01234       theta_default_ = ST::one()/(2*ST::one());
01235       theta_ = theta_default_;
01236       this->setupData();
01237 
01238       RCP<ParameterList> validPL = Teuchos::parameterList();
01239       validPL->set("Description","",this->getMyDescription());
01240       validPL->set<double>("theta",theta_default_,
01241         "Valid values are 0 <= theta <= 1, where theta = 0 "
01242         "implies Forward Euler, theta = 1/2 implies midpoint "
01243         "method, and theta = 1 implies Backward Euler. "
01244         "For theta != 1/2, this method is first-order accurate, "
01245         "and with theta = 1/2, it is second-order accurate.  "
01246         "This method is A-stable, but becomes L-stable with theta=1.");
01247       Teuchos::setupVerboseObjectSublist(&*validPL);
01248       this->setMyValidParameterList(validPL);
01249     }
01250 
01251     void setupData()
01252     {
01253       typedef ScalarTraits<Scalar> ST;
01254       int myNumStages = 1;
01255       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01256       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01257       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01258       myA(0,0) = theta_;
01259       myb(0) = ST::one();
01260       myc(0) = theta_;
01261       this->setMy_A(myA);
01262       this->setMy_b(myb);
01263       this->setMy_c(myc);
01264       if (theta_ == theta_default_) this->setMy_order(2);
01265       else                          this->setMy_order(1);
01266     }
01267 
01268     void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
01269     {
01270       TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
01271       paramList->validateParameters(*this->getValidParameters());
01272       Teuchos::readVerboseObjectSublist(&*paramList,this);
01273       theta_ = paramList->get<double>("theta",theta_default_);
01274       this->setupData();
01275       this->setMyParamList(paramList);
01276     }
01277   private:
01278     Scalar theta_default_;
01279     Scalar theta_;
01280 };
01281 
01282 
01283 template<class Scalar>
01284 class IRK2StageTheta_RKBT :
01285   virtual public RKButcherTableauDefaultBase<Scalar>
01286 {
01287   public:
01288     IRK2StageTheta_RKBT()
01289     {
01290       std::ostringstream myDescription;
01291       myDescription << IRK2StageTheta_name() << "\n"
01292                   << "Computer Methods for ODEs and DAEs\n"
01293                   << "U. M. Ascher and L. R. Petzold\n"
01294                   << "p. 113\n"
01295                   << "c = [  0       1     ]'\n"
01296                   << "A = [  0       0     ]\n"
01297                   << "    [ 1-theta  theta ]\n"
01298                   << "b = [ 1-theta  theta ]'\n"
01299                   << std::endl;
01300 
01301       this->setMyDescription(myDescription.str());
01302       typedef ScalarTraits<Scalar> ST;
01303       theta_default_ = ST::one()/(2*ST::one());
01304       theta_ = theta_default_;
01305       this->setupData();
01306 
01307       RCP<ParameterList> validPL = Teuchos::parameterList();
01308       validPL->set("Description","",this->getMyDescription());
01309       validPL->set<double>("theta",theta_default_,
01310         "Valid values are 0 < theta <= 1, where theta = 0 "
01311         "implies Forward Euler, theta = 1/2 implies trapezoidal "
01312         "method, and theta = 1 implies Backward Euler. "
01313         "For theta != 1/2, this method is first-order accurate, "
01314         "and with theta = 1/2, it is second-order accurate.  "
01315         "This method is A-stable, but becomes L-stable with theta=1.");
01316       Teuchos::setupVerboseObjectSublist(&*validPL);
01317       this->setMyValidParameterList(validPL);
01318     }
01319     void setupData()
01320     {
01321       typedef ScalarTraits<Scalar> ST;
01322       int myNumStages = 2;
01323       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01324       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01325       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01326       Scalar one = ST::one();
01327       Scalar zero = ST::zero();
01328       myA(0,0) = zero;
01329       myA(0,1) = zero;
01330       myA(1,0) = as<Scalar>( one - theta_ );
01331       myA(1,1) = theta_;
01332       myb(0) = as<Scalar>( one - theta_ );
01333       myb(1) = theta_;
01334       myc(0) = theta_;
01335       myc(1) = one;
01336 
01337       this->setMy_A(myA);
01338       this->setMy_b(myb);
01339       this->setMy_c(myc);
01340       TEUCHOS_TEST_FOR_EXCEPTION(
01341         theta_ == zero, std::logic_error,
01342         "'theta' can not be zero, as it makes this IRK stepper explicit.");
01343       if (theta_ == theta_default_) this->setMy_order(2);
01344       else                          this->setMy_order(1);
01345     }
01346     void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
01347     {
01348       TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
01349       paramList->validateParameters(*this->getValidParameters());
01350       Teuchos::readVerboseObjectSublist(&*paramList,this);
01351       theta_ = paramList->get<double>("theta",theta_default_);
01352       this->setupData();
01353       this->setMyParamList(paramList);
01354     }
01355   private:
01356     Scalar theta_default_;
01357     Scalar theta_;
01358 };
01359 
01360 
01361 template<class Scalar>
01362 class Implicit1Stage2ndOrderGauss_RKBT :
01363   virtual public RKButcherTableauDefaultBase<Scalar>
01364 {
01365   public:
01366     Implicit1Stage2ndOrderGauss_RKBT()
01367     {
01368 
01369       std::ostringstream myDescription;
01370       myDescription << Implicit1Stage2ndOrderGauss_name() << "\n"
01371                   << "A-stable\n"
01372                   << "Solving Ordinary Differential Equations II:\n"
01373                   << "Stiff and Differential-Algebraic Problems,\n"
01374                   << "2nd Revised Edition\n"
01375                   << "E. Hairer and G. Wanner\n"
01376                   << "Table 5.2, pg 72\n"
01377                   << "Also:  Implicit midpoint rule\n"
01378                   << "Solving Ordinary Differential Equations I:\n"
01379                   << "Nonstiff Problems, 2nd Revised Edition\n"
01380                   << "E. Hairer, S. P. Norsett, and G. Wanner\n"
01381                   << "Table 7.1, pg 205\n"
01382                   << "c = [ 1/2 ]'\n"
01383                   << "A = [ 1/2 ]\n"
01384                   << "b = [  1  ]'" << std::endl;
01385       typedef ScalarTraits<Scalar> ST;
01386       int myNumStages = 1;
01387       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01388       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01389       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01390       Scalar onehalf = ST::one()/(2*ST::one());
01391       Scalar one = ST::one();
01392       myA(0,0) = onehalf;
01393       myb(0) = one;
01394       myc(0) = onehalf;
01395       this->setMyDescription(myDescription.str());
01396       this->setMy_A(myA);
01397       this->setMy_b(myb);
01398       this->setMy_c(myc);
01399       this->setMy_order(2);
01400     }
01401 };
01402 
01403 
01404 template<class Scalar>
01405 class Implicit2Stage4thOrderGauss_RKBT :
01406   virtual public RKButcherTableauDefaultBase<Scalar>
01407 {
01408   public:
01409     Implicit2Stage4thOrderGauss_RKBT()
01410     {
01411 
01412       std::ostringstream myDescription;
01413       myDescription << Implicit2Stage4thOrderGauss_name() << "\n"
01414                   << "A-stable\n"
01415                   << "Solving Ordinary Differential Equations II:\n"
01416                   << "Stiff and Differential-Algebraic Problems,\n"
01417                   << "2nd Revised Edition\n"
01418                   << "E. Hairer and G. Wanner\n"
01419                   << "Table 5.2, pg 72\n"
01420                   << "c = [ 1/2-sqrt(3)/6  1/2+sqrt(3)/6 ]'\n"
01421                   << "A = [ 1/4            1/4-sqrt(3)/6 ]\n"
01422                   << "    [ 1/4+sqrt(3)/6  1/4           ]\n"
01423                   << "b = [ 1/2            1/2 ]'" << std::endl;
01424       typedef ScalarTraits<Scalar> ST;
01425       int myNumStages = 2;
01426       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01427       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01428       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01429       Scalar one = ST::one();
01430       Scalar onehalf = as<Scalar>(one/(2*one));
01431       Scalar three = as<Scalar>(3*one);
01432       Scalar six = as<Scalar>(6*one);
01433       Scalar onefourth = as<Scalar>(one/(4*one));
01434       Scalar alpha = ST::squareroot(three)/six;
01435 
01436       myA(0,0) = onefourth;
01437       myA(0,1) = onefourth-alpha;
01438       myA(1,0) = onefourth+alpha;
01439       myA(1,1) = onefourth;
01440       myb(0) = onehalf;
01441       myb(1) = onehalf;
01442       myc(0) = onehalf-alpha;
01443       myc(1) = onehalf+alpha;
01444       this->setMyDescription(myDescription.str());
01445       this->setMy_A(myA);
01446       this->setMy_b(myb);
01447       this->setMy_c(myc);
01448       this->setMy_order(4);
01449     }
01450 };
01451 
01452 
01453 template<class Scalar>
01454 class Implicit3Stage6thOrderGauss_RKBT :
01455   virtual public RKButcherTableauDefaultBase<Scalar>
01456 {
01457   public:
01458     Implicit3Stage6thOrderGauss_RKBT()
01459     {
01460 
01461       std::ostringstream myDescription;
01462       myDescription << Implicit3Stage6thOrderGauss_name() << "\n"
01463                   << "A-stable\n"
01464                   << "Solving Ordinary Differential Equations II:\n"
01465                   << "Stiff and Differential-Algebraic Problems,\n"
01466                   << "2nd Revised Edition\n"
01467                   << "E. Hairer and G. Wanner\n"
01468                   << "Table 5.2, pg 72\n"
01469                   << "c = [ 1/2-sqrt(15)/10   1/2              1/2+sqrt(15)/10  ]'\n"
01470                   << "A = [ 5/36              2/9-sqrt(15)/15  5/36-sqrt(15)/30 ]\n"
01471                   << "    [ 5/36+sqrt(15)/24  2/9              5/36-sqrt(15)/24 ]\n"
01472                   << "    [ 5/36+sqrt(15)/30  2/9+sqrt(15)/15  5/36             ]\n"
01473                   << "b = [ 5/18              4/9              5/18             ]'" << std::endl;
01474       typedef ScalarTraits<Scalar> ST;
01475       int myNumStages = 3;
01476       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01477       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01478       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01479       Scalar one = ST::one();
01480       Scalar ten = as<Scalar>(10*one);
01481       Scalar fifteen = as<Scalar>(15*one);
01482       Scalar twentyfour = as<Scalar>(24*one);
01483       Scalar thirty = as<Scalar>(30*one);
01484       Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
01485       Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
01486       Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
01487       Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
01488 
01489       myA(0,0) = as<Scalar>(5*one/(36*one));
01490       myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
01491       myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
01492       myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
01493       myA(1,1) = as<Scalar>(2*one/(9*one));
01494       myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
01495       myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
01496       myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
01497       myA(2,2) = as<Scalar>(5*one/(36*one));
01498       myb(0) = as<Scalar>(5*one/(18*one));
01499       myb(1) = as<Scalar>(4*one/(9*one));
01500       myb(2) = as<Scalar>(5*one/(18*one));
01501       myc(0) = as<Scalar>(one/(2*one))-sqrt15over10;
01502       myc(1) = as<Scalar>(one/(2*one));
01503       myc(2) = as<Scalar>(one/(2*one))+sqrt15over10;
01504       this->setMyDescription(myDescription.str());
01505       this->setMy_A(myA);
01506       this->setMy_b(myb);
01507       this->setMy_c(myc);
01508       this->setMy_order(6);
01509     }
01510 };
01511 
01512 
01513 template<class Scalar>
01514 class Implicit1Stage1stOrderRadauA_RKBT :
01515   virtual public RKButcherTableauDefaultBase<Scalar>
01516 {
01517   public:
01518     Implicit1Stage1stOrderRadauA_RKBT()
01519     {
01520 
01521       std::ostringstream myDescription;
01522       myDescription << Implicit1Stage1stOrderRadauA_name() << "\n"
01523                   << "A-stable\n"
01524                   << "Solving Ordinary Differential Equations II:\n"
01525                   << "Stiff and Differential-Algebraic Problems,\n"
01526                   << "2nd Revised Edition\n"
01527                   << "E. Hairer and G. Wanner\n"
01528                   << "Table 5.3, pg 73\n"
01529                   << "c = [ 0 ]'\n"
01530                   << "A = [ 1 ]\n"
01531                   << "b = [ 1 ]'" << std::endl;
01532       typedef ScalarTraits<Scalar> ST;
01533       int myNumStages = 1;
01534       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01535       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01536       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01537       Scalar one = ST::one();
01538       Scalar zero = ST::zero();
01539       myA(0,0) = one;
01540       myb(0) = one;
01541       myc(0) = zero;
01542       this->setMyDescription(myDescription.str());
01543       this->setMy_A(myA);
01544       this->setMy_b(myb);
01545       this->setMy_c(myc);
01546       this->setMy_order(1);
01547     }
01548 };
01549 
01550 
01551 template<class Scalar>
01552 class Implicit2Stage3rdOrderRadauA_RKBT :
01553   virtual public RKButcherTableauDefaultBase<Scalar>
01554 {
01555   public:
01556     Implicit2Stage3rdOrderRadauA_RKBT()
01557     {
01558 
01559       std::ostringstream myDescription;
01560       myDescription << Implicit2Stage3rdOrderRadauA_name() << "\n"
01561                   << "A-stable\n"
01562                   << "Solving Ordinary Differential Equations II:\n"
01563                   << "Stiff and Differential-Algebraic Problems,\n"
01564                   << "2nd Revised Edition\n"
01565                   << "E. Hairer and G. Wanner\n"
01566                   << "Table 5.3, pg 73\n"
01567                   << "c = [  0    2/3 ]'\n"
01568                   << "A = [ 1/4  -1/4 ]\n"
01569                   << "    [ 1/4  5/12 ]\n"
01570                   << "b = [ 1/4  3/4  ]'" << std::endl;
01571       typedef ScalarTraits<Scalar> ST;
01572       int myNumStages = 2;
01573       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01574       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01575       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01576       Scalar zero = ST::zero();
01577       Scalar one = ST::one();
01578       myA(0,0) = as<Scalar>(one/(4*one));
01579       myA(0,1) = as<Scalar>(-one/(4*one));
01580       myA(1,0) = as<Scalar>(one/(4*one));
01581       myA(1,1) = as<Scalar>(5*one/(12*one));
01582       myb(0) = as<Scalar>(one/(4*one));
01583       myb(1) = as<Scalar>(3*one/(4*one));
01584       myc(0) = zero;
01585       myc(1) = as<Scalar>(2*one/(3*one));
01586       this->setMyDescription(myDescription.str());
01587       this->setMy_A(myA);
01588       this->setMy_b(myb);
01589       this->setMy_c(myc);
01590       this->setMy_order(3);
01591     }
01592 };
01593 
01594 
01595 template<class Scalar>
01596 class Implicit3Stage5thOrderRadauA_RKBT :
01597   virtual public RKButcherTableauDefaultBase<Scalar>
01598 {
01599   public:
01600     Implicit3Stage5thOrderRadauA_RKBT()
01601     {
01602 
01603       std::ostringstream myDescription;
01604       myDescription << Implicit3Stage5thOrderRadauA_name() << "\n"
01605                   << "A-stable\n"
01606                   << "Solving Ordinary Differential Equations II:\n"
01607                   << "Stiff and Differential-Algebraic Problems,\n"
01608                   << "2nd Revised Edition\n"
01609                   << "E. Hairer and G. Wanner\n"
01610                   << "Table 5.4, pg 73\n"
01611                   << "c = [  0   (6-sqrt(6))/10       (6+sqrt(6))/10      ]'\n"
01612                   << "A = [ 1/9  (-1-sqrt(6))/18      (-1+sqrt(6))/18     ]\n"
01613                   << "    [ 1/9  (88+7*sqrt(6))/360   (88-43*sqrt(6))/360 ]\n"
01614                   << "    [ 1/9  (88+43*sqrt(6))/360  (88-7*sqrt(6))/360  ]\n"
01615                   << "b = [ 1/9  (16+sqrt(6))/36      (16-sqrt(6))/36     ]'" << std::endl;
01616       typedef ScalarTraits<Scalar> ST;
01617       int myNumStages = 3;
01618       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01619       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01620       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01621       Scalar zero = ST::zero();
01622       Scalar one = ST::one();
01623       myA(0,0) = as<Scalar>(one/(9*one));
01624       myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
01625       myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
01626       myA(1,0) = as<Scalar>(one/(9*one));
01627       myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
01628       myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
01629       myA(2,0) = as<Scalar>(one/(9*one));
01630       myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
01631       myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
01632       myb(0) = as<Scalar>(one/(9*one));
01633       myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
01634       myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
01635       myc(0) = zero;
01636       myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
01637       myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
01638       this->setMyDescription(myDescription.str());
01639       this->setMy_A(myA);
01640       this->setMy_b(myb);
01641       this->setMy_c(myc);
01642       this->setMy_order(5);
01643     }
01644 };
01645 
01646 
01647 template<class Scalar>
01648 class Implicit1Stage1stOrderRadauB_RKBT :
01649   virtual public RKButcherTableauDefaultBase<Scalar>
01650 {
01651   public:
01652     Implicit1Stage1stOrderRadauB_RKBT()
01653     {
01654 
01655       std::ostringstream myDescription;
01656       myDescription << Implicit1Stage1stOrderRadauB_name() << "\n"
01657                   << "A-stable\n"
01658                   << "Solving Ordinary Differential Equations II:\n"
01659                   << "Stiff and Differential-Algebraic Problems,\n"
01660                   << "2nd Revised Edition\n"
01661                   << "E. Hairer and G. Wanner\n"
01662                   << "Table 5.5, pg 74\n"
01663                   << "c = [ 1 ]'\n"
01664                   << "A = [ 1 ]\n"
01665                   << "b = [ 1 ]'" << std::endl;
01666       typedef ScalarTraits<Scalar> ST;
01667       int myNumStages = 1;
01668       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01669       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01670       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01671       Scalar one = ST::one();
01672       myA(0,0) = one;
01673       myb(0) = one;
01674       myc(0) = one;
01675       this->setMyDescription(myDescription.str());
01676       this->setMy_A(myA);
01677       this->setMy_b(myb);
01678       this->setMy_c(myc);
01679       this->setMy_order(1);
01680     }
01681 };
01682 
01683 
01684 template<class Scalar>
01685 class Implicit2Stage3rdOrderRadauB_RKBT :
01686   virtual public RKButcherTableauDefaultBase<Scalar>
01687 {
01688   public:
01689     Implicit2Stage3rdOrderRadauB_RKBT()
01690     {
01691 
01692       std::ostringstream myDescription;
01693       myDescription << Implicit2Stage3rdOrderRadauB_name() << "\n"
01694                   << "A-stable\n"
01695                   << "Solving Ordinary Differential Equations II:\n"
01696                   << "Stiff and Differential-Algebraic Problems,\n"
01697                   << "2nd Revised Edition\n"
01698                   << "E. Hairer and G. Wanner\n"
01699                   << "Table 5.5, pg 74\n"
01700                   << "c = [ 1/3     1   ]'\n"
01701                   << "A = [ 5/12  -1/12 ]\n"
01702                   << "    [ 3/4    1/4  ]\n"
01703                   << "b = [ 3/4    1/4  ]'" << std::endl;
01704       typedef ScalarTraits<Scalar> ST;
01705       int myNumStages = 2;
01706       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01707       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01708       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01709       Scalar one = ST::one();
01710       myA(0,0) = as<Scalar>( 5*one/(12*one) );
01711       myA(0,1) = as<Scalar>( -one/(12*one) );
01712       myA(1,0) = as<Scalar>( 3*one/(4*one) );
01713       myA(1,1) = as<Scalar>( one/(4*one) );
01714       myb(0) = as<Scalar>( 3*one/(4*one) );
01715       myb(1) = as<Scalar>( one/(4*one) );
01716       myc(0) = as<Scalar>( one/(3*one) );
01717       myc(1) = one;
01718       this->setMyDescription(myDescription.str());
01719       this->setMy_A(myA);
01720       this->setMy_b(myb);
01721       this->setMy_c(myc);
01722       this->setMy_order(3);
01723     }
01724 };
01725 
01726 
01727 template<class Scalar>
01728 class Implicit3Stage5thOrderRadauB_RKBT :
01729   virtual public RKButcherTableauDefaultBase<Scalar>
01730 {
01731   public:
01732     Implicit3Stage5thOrderRadauB_RKBT()
01733     {
01734 
01735       std::ostringstream myDescription;
01736       myDescription << Implicit3Stage5thOrderRadauB_name() << "\n"
01737         << "A-stable\n"
01738         << "Solving Ordinary Differential Equations II:\n"
01739         << "Stiff and Differential-Algebraic Problems,\n"
01740         << "2nd Revised Edition\n"
01741         << "E. Hairer and G. Wanner\n"
01742         << "Table 5.6, pg 74\n"
01743         << "c = [ (4-sqrt(6))/10          (4+sqrt(6))/10          1    ]'\n"
01744         << "A = [ A1 A2 A3 ]\n"
01745         << "      A1 = [ (88-7*sqrt(6))/360     ]\n"
01746         << "           [ (296+169*sqrt(6))/1800 ]\n"
01747         << "           [ (16-sqrt(6))/36        ]\n"
01748         << "      A2 = [ (296-169*sqrt(6))/1800 ]\n"
01749         << "           [ (88+7*sqrt(6))/360     ]\n"
01750         << "           [ (16+sqrt(6))/36        ]\n"
01751         << "      A3 = [ (-2+3*sqrt(6))/225 ]\n"
01752         << "           [ (-2-3*sqrt(6))/225 ]\n"
01753         << "           [ 1/9                ]\n"
01754         << "b = [ (16-sqrt(6))/36         (16+sqrt(6))/36         1/9 ]'"
01755         << std::endl;
01756       typedef ScalarTraits<Scalar> ST;
01757       int myNumStages = 3;
01758       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01759       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01760       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01761       Scalar one = ST::one();
01762       myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
01763       myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
01764       myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
01765       myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
01766       myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
01767       myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
01768       myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
01769       myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
01770       myA(2,2) = as<Scalar>( one/(9*one) );
01771       myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
01772       myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
01773       myb(2) = as<Scalar>( one/(9*one) );
01774       myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
01775       myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
01776       myc(2) = one;
01777       this->setMyDescription(myDescription.str());
01778       this->setMy_A(myA);
01779       this->setMy_b(myb);
01780       this->setMy_c(myc);
01781       this->setMy_order(5);
01782     }
01783 };
01784 
01785 
01786 template<class Scalar>
01787 class Implicit2Stage2ndOrderLobattoA_RKBT :
01788   virtual public RKButcherTableauDefaultBase<Scalar>
01789 {
01790   public:
01791     Implicit2Stage2ndOrderLobattoA_RKBT()
01792     {
01793 
01794       std::ostringstream myDescription;
01795       myDescription << Implicit2Stage2ndOrderLobattoA_name() << "\n"
01796                   << "A-stable\n"
01797                   << "Solving Ordinary Differential Equations II:\n"
01798                   << "Stiff and Differential-Algebraic Problems,\n"
01799                   << "2nd Revised Edition\n"
01800                   << "E. Hairer and G. Wanner\n"
01801                   << "Table 5.7, pg 75\n"
01802                   << "c = [  0    1   ]'\n"
01803                   << "A = [  0    0   ]\n"
01804                   << "    [ 1/2  1/2  ]\n"
01805                   << "b = [ 1/2  1/2  ]'" << std::endl;
01806       typedef ScalarTraits<Scalar> ST;
01807       int myNumStages = 2;
01808       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01809       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01810       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01811       Scalar zero = ST::zero();
01812       Scalar one = ST::one();
01813       myA(0,0) = zero;
01814       myA(0,1) = zero;
01815       myA(1,0) = as<Scalar>( one/(2*one) );
01816       myA(1,1) = as<Scalar>( one/(2*one) );
01817       myb(0) = as<Scalar>( one/(2*one) );
01818       myb(1) = as<Scalar>( one/(2*one) );
01819       myc(0) = zero;
01820       myc(1) = one;
01821       this->setMyDescription(myDescription.str());
01822       this->setMy_A(myA);
01823       this->setMy_b(myb);
01824       this->setMy_c(myc);
01825       this->setMy_order(2);
01826     }
01827 };
01828 
01829 
01830 template<class Scalar>
01831 class Implicit3Stage4thOrderLobattoA_RKBT :
01832   virtual public RKButcherTableauDefaultBase<Scalar>
01833 {
01834   public:
01835     Implicit3Stage4thOrderLobattoA_RKBT()
01836     {
01837 
01838       std::ostringstream myDescription;
01839       myDescription << Implicit3Stage4thOrderLobattoA_name() << "\n"
01840                   << "A-stable\n"
01841                   << "Solving Ordinary Differential Equations II:\n"
01842                   << "Stiff and Differential-Algebraic Problems,\n"
01843                   << "2nd Revised Edition\n"
01844                   << "E. Hairer and G. Wanner\n"
01845                   << "Table 5.7, pg 75\n"
01846                   << "c = [  0    1/2    1  ]'\n"
01847                   << "A = [  0     0     0   ]\n"
01848                   << "    [ 5/24  1/3  -1/24  ]\n"
01849                   << "    [ 1/6   2/3   1/6   ]\n"
01850                   << "b = [ 1/6   2/3   1/6   ]'" << std::endl;
01851       typedef ScalarTraits<Scalar> ST;
01852       int myNumStages = 3;
01853       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01854       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01855       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01856       Scalar zero = ST::zero();
01857       Scalar one = ST::one();
01858       myA(0,0) = zero;
01859       myA(0,1) = zero;
01860       myA(0,2) = zero;
01861       myA(1,0) = as<Scalar>( (5*one)/(24*one) );
01862       myA(1,1) = as<Scalar>( (one)/(3*one) );
01863       myA(1,2) = as<Scalar>( (-one)/(24*one) );
01864       myA(2,0) = as<Scalar>( (one)/(6*one) );
01865       myA(2,1) = as<Scalar>( (2*one)/(3*one) );
01866       myA(2,2) = as<Scalar>( (1*one)/(6*one) );
01867       myb(0) = as<Scalar>( (one)/(6*one) );
01868       myb(1) = as<Scalar>( (2*one)/(3*one) );
01869       myb(2) = as<Scalar>( (1*one)/(6*one) );
01870       myc(0) = zero;
01871       myc(1) = as<Scalar>( one/(2*one) );
01872       myc(2) = one;
01873       this->setMyDescription(myDescription.str());
01874       this->setMy_A(myA);
01875       this->setMy_b(myb);
01876       this->setMy_c(myc);
01877       this->setMy_order(4);
01878     }
01879 };
01880 
01881 
01882 template<class Scalar>
01883 class Implicit4Stage6thOrderLobattoA_RKBT :
01884   virtual public RKButcherTableauDefaultBase<Scalar>
01885 {
01886   public:
01887     Implicit4Stage6thOrderLobattoA_RKBT()
01888     {
01889 
01890       using Teuchos::as;
01891       typedef Teuchos::ScalarTraits<Scalar> ST;
01892 
01893       std::ostringstream myDescription;
01894       myDescription << Implicit4Stage6thOrderLobattoA_name() << "\n"
01895         << "A-stable\n"
01896         << "Solving Ordinary Differential Equations II:\n"
01897         << "Stiff and Differential-Algebraic Problems,\n"
01898         << "2nd Revised Edition\n"
01899         << "E. Hairer and G. Wanner\n"
01900         << "Table 5.8, pg 75\n"
01901         << "c = [ 0  (5-sqrt(5))/10  (5+sqrt(5))/10  1 ]'\n"
01902         << "A = [ A1  A2  A3  A4 ]\n"
01903         << "      A1 = [ 0               ]\n"
01904         << "           [ (11+sqrt(5)/120 ]\n"
01905         << "           [ (11-sqrt(5)/120 ]\n"
01906         << "           [ 1/12            ]\n"
01907         << "      A2 = [ 0                    ]\n"
01908         << "           [ (25-sqrt(5))/120     ]\n"
01909         << "           [ (25+13*sqrt(5))/120  ]\n"
01910         << "           [ 5/12                 ]\n"
01911         << "      A3 = [ 0                   ]\n"
01912         << "           [ (25-13*sqrt(5))/120 ]\n"
01913         << "           [ (25+sqrt(5))/120    ]\n"
01914         << "           [ 5/12                ]\n"
01915         << "      A4 = [ 0                ]\n"
01916         << "           [ (-1+sqrt(5))/120 ]\n"
01917         << "           [ (-1-sqrt(5))/120 ]\n"
01918         << "           [ 1/12             ]\n"
01919         << "b = [ 1/12  5/12  5/12   1/12 ]'" << std::endl;
01920       typedef ScalarTraits<Scalar> ST;
01921       int myNumStages = 4;
01922       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01923       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01924       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01925       Scalar zero = ST::zero();
01926       Scalar one = ST::one();
01927       myA(0,0) = zero;
01928       myA(0,1) = zero;
01929       myA(0,2) = zero;
01930       myA(0,3) = zero;
01931       myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
01932       myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
01933       myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
01934       myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
01935       myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
01936       myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
01937       myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
01938       myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
01939       myA(3,0) = as<Scalar>( one/(12*one) );
01940       myA(3,1) = as<Scalar>( 5*one/(12*one) );
01941       myA(3,2) = as<Scalar>( 5*one/(12*one) );
01942       myA(3,3) = as<Scalar>( one/(12*one) );
01943       myb(0) = as<Scalar>( one/(12*one) );
01944       myb(1) = as<Scalar>( 5*one/(12*one) );
01945       myb(2) = as<Scalar>( 5*one/(12*one) );
01946       myb(3) = as<Scalar>( one/(12*one) );
01947       myc(0) = zero;
01948       myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
01949       myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
01950       myc(3) = one;
01951       this->setMyDescription(myDescription.str());
01952       this->setMy_A(myA);
01953       this->setMy_b(myb);
01954       this->setMy_c(myc);
01955       this->setMy_order(6);
01956     }
01957 };
01958 
01959 
01960 template<class Scalar>
01961 class Implicit2Stage2ndOrderLobattoB_RKBT :
01962   virtual public RKButcherTableauDefaultBase<Scalar>
01963 {
01964   public:
01965     Implicit2Stage2ndOrderLobattoB_RKBT()
01966     {
01967 
01968       std::ostringstream myDescription;
01969       myDescription << Implicit2Stage2ndOrderLobattoB_name() << "\n"
01970                   << "A-stable\n"
01971                   << "Solving Ordinary Differential Equations II:\n"
01972                   << "Stiff and Differential-Algebraic Problems,\n"
01973                   << "2nd Revised Edition\n"
01974                   << "E. Hairer and G. Wanner\n"
01975                   << "Table 5.9, pg 76\n"
01976                   << "c = [  0    1   ]'\n"
01977                   << "A = [ 1/2   0   ]\n"
01978                   << "    [ 1/2   0   ]\n"
01979                   << "b = [ 1/2  1/2  ]'" << std::endl;
01980       typedef ScalarTraits<Scalar> ST;
01981       int myNumStages = 2;
01982       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
01983       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
01984       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
01985       Scalar zero = ST::zero();
01986       Scalar one = ST::one();
01987       myA(0,0) = as<Scalar>( one/(2*one) );
01988       myA(0,1) = zero;
01989       myA(1,0) = as<Scalar>( one/(2*one) );
01990       myA(1,1) = zero;
01991       myb(0) = as<Scalar>( one/(2*one) );
01992       myb(1) = as<Scalar>( one/(2*one) );
01993       myc(0) = zero;
01994       myc(1) = one;
01995       this->setMyDescription(myDescription.str());
01996       this->setMy_A(myA);
01997       this->setMy_b(myb);
01998       this->setMy_c(myc);
01999       this->setMy_order(2);
02000     }
02001 };
02002 
02003 
02004 template<class Scalar>
02005 class Implicit3Stage4thOrderLobattoB_RKBT :
02006   virtual public RKButcherTableauDefaultBase<Scalar>
02007 {
02008   public:
02009     Implicit3Stage4thOrderLobattoB_RKBT()
02010     {
02011 
02012       std::ostringstream myDescription;
02013       myDescription << Implicit3Stage4thOrderLobattoB_name() << "\n"
02014                   << "A-stable\n"
02015                   << "Solving Ordinary Differential Equations II:\n"
02016                   << "Stiff and Differential-Algebraic Problems,\n"
02017                   << "2nd Revised Edition\n"
02018                   << "E. Hairer and G. Wanner\n"
02019                   << "Table 5.9, pg 76\n"
02020                   << "c = [  0    1/2    1   ]'\n"
02021                   << "A = [ 1/6  -1/6    0   ]\n"
02022                   << "    [ 1/6   1/3    0   ]\n"
02023                   << "    [ 1/6   5/6    0   ]\n"
02024                   << "b = [ 1/6   2/3   1/6  ]'" << std::endl;
02025       typedef ScalarTraits<Scalar> ST;
02026       int myNumStages = 3;
02027       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02028       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02029       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02030       Scalar zero = ST::zero();
02031       Scalar one = ST::one();
02032       myA(0,0) = as<Scalar>( one/(6*one) );
02033       myA(0,1) = as<Scalar>( -one/(6*one) );
02034       myA(0,2) = zero;
02035       myA(1,0) = as<Scalar>( one/(6*one) );
02036       myA(1,1) = as<Scalar>( one/(3*one) );
02037       myA(1,2) = zero;
02038       myA(2,0) = as<Scalar>( one/(6*one) );
02039       myA(2,1) = as<Scalar>( 5*one/(6*one) );
02040       myA(2,2) = zero;
02041       myb(0) = as<Scalar>( one/(6*one) );
02042       myb(1) = as<Scalar>( 2*one/(3*one) );
02043       myb(2) = as<Scalar>( one/(6*one) );
02044       myc(0) = zero;
02045       myc(1) = as<Scalar>( one/(2*one) );
02046       myc(2) = one;
02047       this->setMyDescription(myDescription.str());
02048       this->setMy_A(myA);
02049       this->setMy_b(myb);
02050       this->setMy_c(myc);
02051       this->setMy_order(4);
02052     }
02053 };
02054 
02055 
02056 template<class Scalar>
02057 class Implicit4Stage6thOrderLobattoB_RKBT :
02058   virtual public RKButcherTableauDefaultBase<Scalar>
02059 {
02060   public:
02061     Implicit4Stage6thOrderLobattoB_RKBT()
02062     {
02063 
02064       std::ostringstream myDescription;
02065       myDescription << Implicit4Stage6thOrderLobattoB_name() << "\n"
02066                   << "A-stable\n"
02067                   << "Solving Ordinary Differential Equations II:\n"
02068                   << "Stiff and Differential-Algebraic Problems,\n"
02069                   << "2nd Revised Edition\n"
02070                   << "E. Hairer and G. Wanner\n"
02071                   << "Table 5.10, pg 76\n"
02072                   << "c = [ 0     (5-sqrt(5))/10       (5+sqrt(5))/10       1     ]'\n"
02073                   << "A = [ 1/12  (-1-sqrt(5))/24      (-1+sqrt(5))/24      0     ]\n"
02074                   << "    [ 1/12  (25+sqrt(5))/120     (25-13*sqrt(5))/120  0     ]\n"
02075                   << "    [ 1/12  (25+13*sqrt(5))/120  (25-sqrt(5))/120     0     ]\n"
02076                   << "    [ 1/12  (11-sqrt(5))/24      (11+sqrt(5))/24      0     ]\n"
02077                   << "b = [ 1/12  5/12                 5/12                 1/12  ]'" << std::endl;
02078       typedef ScalarTraits<Scalar> ST;
02079       int myNumStages = 4;
02080       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02081       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02082       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02083       Scalar zero = ST::zero();
02084       Scalar one = ST::one();
02085       myA(0,0) = as<Scalar>( one/(12*one) );
02086       myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
02087       myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
02088       myA(0,3) = zero;
02089       myA(1,0) = as<Scalar>( one/(12*one) );
02090       myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
02091       myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
02092       myA(1,3) = zero;
02093       myA(2,0) = as<Scalar>( one/(12*one) );
02094       myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
02095       myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
02096       myA(2,3) = zero;
02097       myA(3,0) = as<Scalar>( one/(12*one) );
02098       myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
02099       myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
02100       myA(3,3) = zero;
02101       myb(0) = as<Scalar>( one/(12*one) );
02102       myb(1) = as<Scalar>( 5*one/(12*one) );
02103       myb(2) = as<Scalar>( 5*one/(12*one) );
02104       myb(3) = as<Scalar>( one/(12*one) );
02105       myc(0) = zero;
02106       myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
02107       myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
02108       myc(3) = one;
02109       this->setMyDescription(myDescription.str());
02110       this->setMy_A(myA);
02111       this->setMy_b(myb);
02112       this->setMy_c(myc);
02113       this->setMy_order(6);
02114     }
02115 };
02116 
02117 
02118 template<class Scalar>
02119 class Implicit2Stage2ndOrderLobattoC_RKBT :
02120   virtual public RKButcherTableauDefaultBase<Scalar>
02121 {
02122   public:
02123     Implicit2Stage2ndOrderLobattoC_RKBT()
02124     {
02125 
02126       std::ostringstream myDescription;
02127       myDescription << Implicit2Stage2ndOrderLobattoC_name() << "\n"
02128                   << "A-stable\n"
02129                   << "Solving Ordinary Differential Equations II:\n"
02130                   << "Stiff and Differential-Algebraic Problems,\n"
02131                   << "2nd Revised Edition\n"
02132                   << "E. Hairer and G. Wanner\n"
02133                   << "Table 5.11, pg 76\n"
02134                   << "c = [  0    1   ]'\n"
02135                   << "A = [ 1/2 -1/2  ]\n"
02136                   << "    [ 1/2  1/2  ]\n"
02137                   << "b = [ 1/2  1/2  ]'" << std::endl;
02138       typedef ScalarTraits<Scalar> ST;
02139       int myNumStages = 2;
02140       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02141       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02142       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02143       Scalar zero = ST::zero();
02144       Scalar one = ST::one();
02145       myA(0,0) = as<Scalar>( one/(2*one) );
02146       myA(0,1) = as<Scalar>( -one/(2*one) );
02147       myA(1,0) = as<Scalar>( one/(2*one) );
02148       myA(1,1) = as<Scalar>( one/(2*one) );
02149       myb(0) = as<Scalar>( one/(2*one) );
02150       myb(1) = as<Scalar>( one/(2*one) );
02151       myc(0) = zero;
02152       myc(1) = one;
02153       this->setMyDescription(myDescription.str());
02154       this->setMy_A(myA);
02155       this->setMy_b(myb);
02156       this->setMy_c(myc);
02157       this->setMy_order(2);
02158     }
02159 };
02160 
02161 
02162 template<class Scalar>
02163 class Implicit3Stage4thOrderLobattoC_RKBT :
02164   virtual public RKButcherTableauDefaultBase<Scalar>
02165 {
02166   public:
02167     Implicit3Stage4thOrderLobattoC_RKBT()
02168     {
02169 
02170       std::ostringstream myDescription;
02171       myDescription << Implicit3Stage4thOrderLobattoC_name() << "\n"
02172                   << "A-stable\n"
02173                   << "Solving Ordinary Differential Equations II:\n"
02174                   << "Stiff and Differential-Algebraic Problems,\n"
02175                   << "2nd Revised Edition\n"
02176                   << "E. Hairer and G. Wanner\n"
02177                   << "Table 5.11, pg 76\n"
02178                   << "c = [  0    1/2    1   ]'\n"
02179                   << "A = [ 1/6  -1/3   1/6  ]\n"
02180                   << "    [ 1/6   5/12 -1/12 ]\n"
02181                   << "    [ 1/6   2/3   1/6  ]\n"
02182                   << "b = [ 1/6   2/3   1/6  ]'" << std::endl;
02183       typedef ScalarTraits<Scalar> ST;
02184       int myNumStages = 3;
02185       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02186       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02187       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02188       Scalar zero = ST::zero();
02189       Scalar one = ST::one();
02190       myA(0,0) = as<Scalar>( one/(6*one) );
02191       myA(0,1) = as<Scalar>( -one/(3*one) );
02192       myA(0,2) = as<Scalar>( one/(6*one) );
02193       myA(1,0) = as<Scalar>( one/(6*one) );
02194       myA(1,1) = as<Scalar>( 5*one/(12*one) );
02195       myA(1,2) = as<Scalar>( -one/(12*one) );
02196       myA(2,0) = as<Scalar>( one/(6*one) );
02197       myA(2,1) = as<Scalar>( 2*one/(3*one) );
02198       myA(2,2) = as<Scalar>( one/(6*one) );
02199       myb(0) = as<Scalar>( one/(6*one) );
02200       myb(1) = as<Scalar>( 2*one/(3*one) );
02201       myb(2) = as<Scalar>( one/(6*one) );
02202       myc(0) = zero;
02203       myc(1) = as<Scalar>( one/(2*one) );
02204       myc(2) = one;
02205       this->setMyDescription(myDescription.str());
02206       this->setMy_A(myA);
02207       this->setMy_b(myb);
02208       this->setMy_c(myc);
02209       this->setMy_order(4);
02210     }
02211 };
02212 
02213 
02214 template<class Scalar>
02215 class Implicit4Stage6thOrderLobattoC_RKBT :
02216   virtual public RKButcherTableauDefaultBase<Scalar>
02217 {
02218   public:
02219     Implicit4Stage6thOrderLobattoC_RKBT()
02220     {
02221 
02222       std::ostringstream myDescription;
02223       myDescription << Implicit4Stage6thOrderLobattoC_name() << "\n"
02224                   << "A-stable\n"
02225                   << "Solving Ordinary Differential Equations II:\n"
02226                   << "Stiff and Differential-Algebraic Problems,\n"
02227                   << "2nd Revised Edition\n"
02228                   << "E. Hairer and G. Wanner\n"
02229                   << "Table 5.12, pg 76\n"
02230                   << "c = [ 0     (5-sqrt(5))/10       (5+sqrt(5))/10       1          ]'\n"
02231                   << "A = [ 1/12  -sqrt(5)/12          sqrt(5)/12          -1/12       ]\n"
02232                   << "    [ 1/12  1/4                  (10-7*sqrt(5))/60   sqrt(5)/60  ]\n"
02233                   << "    [ 1/12  (10+7*sqrt(5))/60    1/4                 -sqrt(5)/60 ]\n"
02234                   << "    [ 1/12  5/12                 5/12                 1/12       ]\n"
02235                   << "b = [ 1/12  5/12                 5/12                 1/12       ]'" << std::endl;
02236       typedef ScalarTraits<Scalar> ST;
02237       int myNumStages = 4;
02238       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02239       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02240       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02241       Scalar zero = ST::zero();
02242       Scalar one = ST::one();
02243       myA(0,0) = as<Scalar>( one/(12*one) );
02244       myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
02245       myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
02246       myA(0,3) = as<Scalar>( -one/(12*one) );
02247       myA(1,0) = as<Scalar>( one/(12*one) );
02248       myA(1,1) = as<Scalar>( one/(4*one) );
02249       myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
02250       myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
02251       myA(2,0) = as<Scalar>( one/(12*one) );
02252       myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
02253       myA(2,2) = as<Scalar>( one/(4*one) );
02254       myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
02255       myA(3,0) = as<Scalar>( one/(12*one) );
02256       myA(3,1) = as<Scalar>( 5*one/(12*one) );
02257       myA(3,2) = as<Scalar>( 5*one/(12*one) );
02258       myA(3,3) = as<Scalar>( one/(12*one) );
02259       myb(0) = as<Scalar>( one/(12*one) );
02260       myb(1) = as<Scalar>( 5*one/(12*one) );
02261       myb(2) = as<Scalar>( 5*one/(12*one) );
02262       myb(3) = as<Scalar>( one/(12*one) );
02263       myc(0) = zero;
02264       myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
02265       myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
02266       myc(3) = one;
02267       this->setMyDescription(myDescription.str());
02268       this->setMy_A(myA);
02269       this->setMy_b(myb);
02270       this->setMy_c(myc);
02271       this->setMy_order(6);
02272     }
02273 };
02274 
02275 
02276 
02277 template<class Scalar>
02278 class SDIRK5Stage5thOrder_RKBT :
02279   virtual public RKButcherTableauDefaultBase<Scalar>
02280 {
02281   public:
02282     SDIRK5Stage5thOrder_RKBT()
02283     {
02284 
02285       std::ostringstream myDescription;
02286       myDescription << SDIRK5Stage5thOrder_name() << "\n"
02287         << "A-stable\n"
02288         << "Solving Ordinary Differential Equations II:\n"
02289         << "Stiff and Differential-Algebraic Problems,\n"
02290         << "2nd Revised Edition\n"
02291         << "E. Hairer and G. Wanner\n"
02292         << "pg101 \n"
02293         << "c = [ (6-sqrt(6))/10   ]\n"
02294         << "    [ (6+9*sqrt(6))/35 ]\n"
02295         << "    [ 1                ]\n"
02296         << "    [ (4-sqrt(6))/10   ]\n"
02297         << "    [ (4+sqrt(6))/10   ]\n"
02298         << "A = [ A1 A2 A3 A4 A5 ]\n"
02299         << "      A1 = [ (6-sqrt(6))/10               ]\n"
02300         << "           [ (-6+5*sqrt(6))/14            ]\n"
02301         << "           [ (888+607*sqrt(6))/2850       ]\n"
02302         << "           [ (3153-3082*sqrt(6))/14250    ]\n"
02303         << "           [ (-32583+14638*sqrt(6))/71250 ]\n"
02304         << "      A2 = [ 0                           ]\n"
02305         << "           [ (6-sqrt(6))/10              ]\n"
02306         << "           [ (126-161*sqrt(6))/1425      ]\n"
02307         << "           [ (3213+1148*sqrt(6))/28500   ]\n"
02308         << "           [ (-17199+364*sqrt(6))/142500 ]\n"
02309         << "      A3 = [ 0                       ]\n"
02310         << "           [ 0                       ]\n"
02311         << "           [ (6-sqrt(6))/10          ]\n"
02312         << "           [ (-267+88*sqrt(6))/500   ]\n"
02313         << "           [ (1329-544*sqrt(6))/2500 ]\n"
02314         << "      A4 = [ 0                     ]\n"
02315         << "           [ 0                     ]\n"
02316         << "           [ 0                     ]\n"
02317         << "           [ (6-sqrt(6))/10        ]\n"
02318         << "           [ (-96+131*sqrt(6))/625 ]\n"
02319         << "      A5 = [ 0              ]\n"
02320         << "           [ 0              ]\n"
02321         << "           [ 0              ]\n"
02322         << "           [ 0              ]\n"
02323         << "           [ (6-sqrt(6))/10 ]\n"
02324         << "b = [               0 ]\n"
02325         << "    [               0 ]\n"
02326         << "    [             1/9 ]\n"
02327         << "    [ (16-sqrt(6))/36 ]\n"
02328         << "    [ (16+sqrt(6))/36 ]" << std::endl;
02329       typedef ScalarTraits<Scalar> ST;
02330       int myNumStages = 5;
02331       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02332       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02333       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02334       Scalar zero = ST::zero();
02335       Scalar one = ST::one();
02336       Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
02337       Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
02338       myA(0,0) = gamma;
02339       myA(0,1) = zero;
02340       myA(0,2) = zero;
02341       myA(0,3) = zero;
02342       myA(0,4) = zero;
02343 
02344       myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
02345       myA(1,1) = gamma;
02346       myA(1,2) = zero;
02347       myA(1,3) = zero;
02348       myA(1,4) = zero;
02349 
02350       myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
02351       myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
02352       myA(2,2) = gamma;
02353       myA(2,3) = zero;
02354       myA(2,4) = zero;
02355 
02356       myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
02357       myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
02358       myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
02359       myA(3,3) = gamma;
02360       myA(3,4) = zero;
02361 
02362       myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
02363       myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
02364       myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
02365       myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
02366       myA(4,4) = gamma;
02367 
02368       myb(0) = zero;
02369       myb(1) = zero;
02370       myb(2) = as<Scalar>( one/(9*one) );
02371       myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
02372       myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
02373 
02374       myc(0) = gamma;
02375       myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
02376       myc(2) = one;
02377       myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
02378       myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
02379 
02380       this->setMyDescription(myDescription.str());
02381       this->setMy_A(myA);
02382       this->setMy_b(myb);
02383       this->setMy_c(myc);
02384       this->setMy_order(5);
02385     }
02386 };
02387 
02388 
02389 template<class Scalar>
02390 class SDIRK5Stage4thOrder_RKBT :
02391   virtual public RKButcherTableauDefaultBase<Scalar>
02392 {
02393   public:
02394     SDIRK5Stage4thOrder_RKBT()
02395     {
02396 
02397       std::ostringstream myDescription;
02398       myDescription << SDIRK5Stage4thOrder_name() << "\n"
02399         << "L-stable\n"
02400         << "Solving Ordinary Differential Equations II:\n"
02401         << "Stiff and Differential-Algebraic Problems,\n"
02402         << "2nd Revised Edition\n"
02403         << "E. Hairer and G. Wanner\n"
02404         << "pg100 \n"
02405         << "c  = [ 1/4       3/4        11/20   1/2     1   ]'\n"
02406         << "A  = [ 1/4                                      ]\n"
02407         << "     [ 1/2       1/4                            ]\n"
02408         << "     [ 17/50     -1/25      1/4                 ]\n"
02409         << "     [ 371/1360  -137/2720  15/544  1/4         ]\n"
02410         << "     [ 25/24     -49/48     125/16  -85/12  1/4 ]\n"
02411         << "b  = [ 25/24     -49/48     125/16  -85/12  1/4 ]'\n"
02412         << "b' = [ 59/48     -17/96     225/32  -85/12  0   ]'" << std::endl;
02413       typedef ScalarTraits<Scalar> ST;
02414       int myNumStages = 5;
02415       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02416       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02417       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02418       Scalar zero = ST::zero();
02419       Scalar one = ST::one();
02420       Scalar onequarter = as<Scalar>( one/(4*one) );
02421       myA(0,0) = onequarter;
02422       myA(0,1) = zero;
02423       myA(0,2) = zero;
02424       myA(0,3) = zero;
02425       myA(0,4) = zero;
02426 
02427       myA(1,0) = as<Scalar>( one / (2*one) );
02428       myA(1,1) = onequarter;
02429       myA(1,2) = zero;
02430       myA(1,3) = zero;
02431       myA(1,4) = zero;
02432 
02433       myA(2,0) = as<Scalar>( 17*one/(50*one) );
02434       myA(2,1) = as<Scalar>( -one/(25*one) );
02435       myA(2,2) = onequarter;
02436       myA(2,3) = zero;
02437       myA(2,4) = zero;
02438 
02439       myA(3,0) = as<Scalar>( 371*one/(1360*one) );
02440       myA(3,1) = as<Scalar>( -137*one/(2720*one) );
02441       myA(3,2) = as<Scalar>( 15*one/(544*one) );
02442       myA(3,3) = onequarter;
02443       myA(3,4) = zero;
02444 
02445       myA(4,0) = as<Scalar>( 25*one/(24*one) );
02446       myA(4,1) = as<Scalar>( -49*one/(48*one) );
02447       myA(4,2) = as<Scalar>( 125*one/(16*one) );
02448       myA(4,3) = as<Scalar>( -85*one/(12*one) );
02449       myA(4,4) = onequarter;
02450 
02451       myb(0) = as<Scalar>( 25*one/(24*one) );
02452       myb(1) = as<Scalar>( -49*one/(48*one) );
02453       myb(2) = as<Scalar>( 125*one/(16*one) );
02454       myb(3) = as<Scalar>( -85*one/(12*one) );
02455       myb(4) = onequarter;
02456 
02457       /*
02458       // Alternate version
02459       myb(0) = as<Scalar>( 59*one/(48*one) );
02460       myb(1) = as<Scalar>( -17*one/(96*one) );
02461       myb(2) = as<Scalar>( 225*one/(32*one) );
02462       myb(3) = as<Scalar>( -85*one/(12*one) );
02463       myb(4) = zero;
02464       */
02465       myc(0) = onequarter;
02466       myc(1) = as<Scalar>( 3*one/(4*one) );
02467       myc(2) = as<Scalar>( 11*one/(20*one) );
02468       myc(3) = as<Scalar>( one/(2*one) );
02469       myc(4) = one;
02470 
02471       this->setMyDescription(myDescription.str());
02472       this->setMy_A(myA);
02473       this->setMy_b(myb);
02474       this->setMy_c(myc);
02475       this->setMy_order(4);
02476     }
02477 };
02478 
02479 
02480 template<class Scalar>
02481 class SDIRK3Stage4thOrder_RKBT :
02482   virtual public RKButcherTableauDefaultBase<Scalar>
02483 {
02484   public:
02485     SDIRK3Stage4thOrder_RKBT()
02486     {
02487 
02488       std::ostringstream myDescription;
02489       myDescription << SDIRK3Stage4thOrder_name() << "\n"
02490                   << "A-stable\n"
02491                   << "Solving Ordinary Differential Equations II:\n"
02492                   << "Stiff and Differential-Algebraic Problems,\n"
02493                   << "2nd Revised Edition\n"
02494                   << "E. Hairer and G. Wanner\n"
02495                   << "pg100 \n"
02496                   << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
02497                   << "delta = 1/(6*(2*gamma-1)^2)\n"
02498                   << "c = [ gamma      1/2        1-gamma ]'\n"
02499                   << "A = [ gamma                         ]\n"
02500                   << "    [ 1/2-gamma  gamma              ]\n"
02501                   << "    [ 2*gamma    1-4*gamma  gamma   ]\n"
02502                   << "b = [ delta      1-2*delta  delta   ]'" << std::endl;
02503       typedef ScalarTraits<Scalar> ST;
02504       int myNumStages = 3;
02505       Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
02506       Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
02507       Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
02508       Scalar zero = ST::zero();
02509       Scalar one = ST::one();
02510       Scalar pi = as<Scalar>(4*one)*std::atan(one);
02511       Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
02512       Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
02513       myA(0,0) = gamma;
02514       myA(0,1) = zero;
02515       myA(0,2) = zero;
02516 
02517       myA(1,0) = as<Scalar>( one/(2*one) - gamma );
02518       myA(1,1) = gamma;
02519       myA(1,2) = zero;
02520 
02521       myA(2,0) = as<Scalar>( 2*gamma );
02522       myA(2,1) = as<Scalar>( one - 4*gamma );
02523       myA(2,2) = gamma;
02524 
02525       myb(0) = delta;
02526       myb(1) = as<Scalar>( one-2*delta );
02527       myb(2) = delta;
02528 
02529       myc(0) = gamma;
02530       myc(1) = as<Scalar>( one/(2*one) );
02531       myc(2) = as<Scalar>( one - gamma );
02532 
02533       this->setMyDescription(myDescription.str());
02534       this->setMy_A(myA);
02535       this->setMy_b(myb);
02536       this->setMy_c(myc);
02537       this->setMy_order(4);
02538     }
02539 };
02540 
02541 
02542 } // namespace Rythmos
02543 
02544 
02545 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP
 All Classes Functions Variables Typedefs Friends