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