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

Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7