|
Rythmos - Transient Integration for Differential Equations
Version of the Day
|
|
Rythmos - Transient Integration for Differential Equations
Version of the Day
|
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
1.7.6.1
