|
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 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>& new_A) { A_ = new_A; } 00192 void setMy_b(const Teuchos::SerialDenseVector<int,Scalar>& new_b) { b_ = new_b; } 00193 void setMy_c(const Teuchos::SerialDenseVector<int,Scalar>& new_c) { c_ = new_c; } 00194 void setMy_order(const int& new_order) { order_ = new_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 myDescription; 00240 myDescription << 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> myA(1,1); 00246 myA(0,0) = ST::one(); 00247 Teuchos::SerialDenseVector<int,Scalar> myb(1); 00248 myb(0) = ST::one(); 00249 Teuchos::SerialDenseVector<int,Scalar> myc(1); 00250 myc(0) = ST::one(); 00251 00252 this->setMyDescription(myDescription.str()); 00253 this->setMy_A(myA); 00254 this->setMy_b(myb); 00255 this->setMy_c(myc); 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 myDescription; 00271 myDescription << 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> myA(1,1); 00277 Teuchos::SerialDenseVector<int,Scalar> myb(1); 00278 myb(0) = ST::one(); 00279 Teuchos::SerialDenseVector<int,Scalar> myc(1); 00280 00281 this->setMyDescription(myDescription.str()); 00282 this->setMy_A(myA); 00283 this->setMy_b(myb); 00284 this->setMy_c(myc); 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 myDescription; 00298 myDescription << 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 myNumStages = 4; 00317 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00318 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00319 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00320 00321 // Fill A: 00322 myA(0,0) = zero; 00323 myA(0,1) = zero; 00324 myA(0,2) = zero; 00325 myA(0,3) = zero; 00326 00327 myA(1,0) = onehalf; 00328 myA(1,1) = zero; 00329 myA(1,2) = zero; 00330 myA(1,3) = zero; 00331 00332 myA(2,0) = zero; 00333 myA(2,1) = onehalf; 00334 myA(2,2) = zero; 00335 myA(2,3) = zero; 00336 00337 myA(3,0) = zero; 00338 myA(3,1) = zero; 00339 myA(3,2) = one; 00340 myA(3,3) = zero; 00341 00342 // Fill myb: 00343 myb(0) = onesixth; 00344 myb(1) = onethird; 00345 myb(2) = onethird; 00346 myb(3) = onesixth; 00347 00348 // fill b_c_ 00349 myc(0) = zero; 00350 myc(1) = onehalf; 00351 myc(2) = onehalf; 00352 myc(3) = one; 00353 00354 this->setMyDescription(myDescription.str()); 00355 this->setMy_A(myA); 00356 this->setMy_b(myb); 00357 this->setMy_c(myc); 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 myDescription; 00372 myDescription << 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 myNumStages = 4; 00384 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00385 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00386 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA: 00396 myA(0,0) = zero; 00397 myA(0,1) = zero; 00398 myA(0,2) = zero; 00399 myA(0,3) = zero; 00400 00401 myA(1,0) = one_third; 00402 myA(1,1) = zero; 00403 myA(1,2) = zero; 00404 myA(1,3) = zero; 00405 00406 myA(2,0) = as<Scalar>(-one_third); 00407 myA(2,1) = one; 00408 myA(2,2) = zero; 00409 myA(2,3) = zero; 00410 00411 myA(3,0) = one; 00412 myA(3,1) = as<Scalar>(-one); 00413 myA(3,2) = one; 00414 myA(3,3) = zero; 00415 00416 // Fill myb: 00417 myb(0) = one_eighth; 00418 myb(1) = three_eighth; 00419 myb(2) = three_eighth; 00420 myb(3) = one_eighth; 00421 00422 // Fill myc: 00423 myc(0) = zero; 00424 myc(1) = one_third; 00425 myc(2) = two_third; 00426 myc(3) = one; 00427 00428 this->setMyDescription(myDescription.str()); 00429 this->setMy_A(myA); 00430 this->setMy_b(myb); 00431 this->setMy_c(myc); 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 myDescription; 00446 myDescription << 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 myNumStages = 4; 00458 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00459 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00460 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = zero; 00470 myA(0,1) = zero; 00471 myA(0,2) = zero; 00472 myA(0,3) = zero; 00473 00474 myA(1,0) = onehalf; 00475 myA(1,1) = zero; 00476 myA(1,2) = zero; 00477 myA(1,3) = zero; 00478 00479 myA(2,0) = zero; 00480 myA(2,1) = one; 00481 myA(2,2) = zero; 00482 myA(2,3) = zero; 00483 00484 myA(3,0) = zero; 00485 myA(3,1) = zero; 00486 myA(3,2) = one; 00487 myA(3,3) = zero; 00488 00489 // Fill b: 00490 myb(0) = onesixth; 00491 myb(1) = twothirds; 00492 myb(2) = zero; 00493 myb(3) = onesixth; 00494 00495 // Fill myc: 00496 myc(0) = zero; 00497 myc(1) = onehalf; 00498 myc(2) = one; 00499 myc(3) = one; 00500 00501 this->setMyDescription(myDescription.str()); 00502 this->setMy_A(myA); 00503 this->setMy_b(myb); 00504 this->setMy_c(myc); 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 myDescription; 00519 myDescription << 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 myNumStages = 3; 00534 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00535 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00536 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00537 00538 // Fill myA: 00539 myA(0,0) = zero; 00540 myA(0,1) = zero; 00541 myA(0,2) = zero; 00542 00543 myA(1,0) = onehalf; 00544 myA(1,1) = zero; 00545 myA(1,2) = zero; 00546 00547 myA(2,0) = -one; 00548 myA(2,1) = two; 00549 myA(2,2) = zero; 00550 00551 // Fill myb: 00552 myb(0) = onesixth; 00553 myb(1) = foursixth; 00554 myb(2) = onesixth; 00555 00556 // fill b_c_ 00557 myc(0) = zero; 00558 myc(1) = onehalf; 00559 myc(2) = one; 00560 00561 this->setMyDescription(myDescription.str()); 00562 this->setMy_A(myA); 00563 this->setMy_b(myb); 00564 this->setMy_c(myc); 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 myDescription; 00578 myDescription << 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 myNumStages = 3; 00596 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00597 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00598 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00599 00600 // Fill myA: 00601 myA(0,0) = zero; 00602 myA(0,1) = zero; 00603 myA(0,2) = zero; 00604 00605 myA(1,0) = onethird; 00606 myA(1,1) = zero; 00607 myA(1,2) = zero; 00608 00609 myA(2,0) = zero; 00610 myA(2,1) = twothirds; 00611 myA(2,2) = zero; 00612 00613 // Fill myb: 00614 myb(0) = onefourth; 00615 myb(1) = zero; 00616 myb(2) = threefourths; 00617 00618 // fill b_c_ 00619 myc(0) = zero; 00620 myc(1) = onethird; 00621 myc(2) = twothirds; 00622 00623 this->setMyDescription(myDescription.str()); 00624 this->setMy_A(myA); 00625 this->setMy_b(myb); 00626 this->setMy_c(myc); 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 myDescription; 00640 myDescription << 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 myNumStages = 2; 00654 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00655 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00656 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00657 00658 // Fill myA: 00659 myA(0,0) = zero; 00660 myA(0,1) = zero; 00661 00662 myA(1,0) = onehalf; 00663 myA(1,1) = zero; 00664 00665 // Fill myb: 00666 myb(0) = zero; 00667 myb(1) = one; 00668 00669 // fill b_c_ 00670 myc(0) = zero; 00671 myc(1) = onehalf; 00672 00673 this->setMyDescription(myDescription.str()); 00674 this->setMy_A(myA); 00675 this->setMy_b(myb); 00676 this->setMy_c(myc); 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 myDescription; 00692 myDescription << 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(myDescription.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 myNumStages = 2; 00717 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00718 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00719 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = gamma; 00724 myA(0,1) = zero; 00725 myA(1,0) = as<Scalar>( one - 2*gamma ); 00726 myA(1,1) = gamma; 00727 myb(0) = as<Scalar>( one/(2*one) ); 00728 myb(1) = as<Scalar>( one/(2*one) ); 00729 myc(0) = gamma; 00730 myc(1) = as<Scalar>( one - gamma ); 00731 00732 this->setMy_A(myA); 00733 this->setMy_b(myb); 00734 this->setMy_c(myc); 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 myDescription; 00761 myDescription << 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 myNumStages = 2; 00772 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00773 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00774 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00775 Scalar one = ST::one(); 00776 Scalar zero = ST::zero(); 00777 myA(0,0) = zero; 00778 myA(0,1) = zero; 00779 myA(1,0) = as<Scalar>( one/(3*one) ); 00780 myA(1,1) = as<Scalar>( one/(3*one) ); 00781 myb(0) = as<Scalar>( one/(4*one) ); 00782 myb(1) = as<Scalar>( 3*one/(4*one) ); 00783 myc(0) = zero; 00784 myc(1) = as<Scalar>( 2*one/(3*one) ); 00785 this->setMyDescription(myDescription.str()); 00786 this->setMy_A(myA); 00787 this->setMy_b(myb); 00788 this->setMy_c(myc); 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 myDescription; 00803 myDescription << 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 myNumStages = 3; 00815 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00816 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00817 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00818 Scalar one = ST::one(); 00819 myA(0,0) = as<Scalar>( 5*one/(36*one) ); 00820 myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) ); 00821 myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) ); 00822 myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) ); 00823 myA(1,1) = as<Scalar>( 2*one/(9*one) ); 00824 myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) ); 00825 myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) ); 00826 myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) ); 00827 myA(2,2) = as<Scalar>( 5*one/(36*one) ); 00828 myb(0) = as<Scalar>( 5*one/(18*one) ); 00829 myb(1) = as<Scalar>( 4*one/(9*one) ); 00830 myb(2) = as<Scalar>( 5*one/(18*one) ); 00831 myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) ); 00832 myc(1) = as<Scalar>( one/(2*one) ); 00833 myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) ); 00834 this->setMyDescription(myDescription.str()); 00835 this->setMy_A(myA); 00836 this->setMy_b(myb); 00837 this->setMy_c(myc); 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 myDescription; 00852 myDescription << 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 myNumStages = 4; 00875 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00876 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00877 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = w1; 00891 myA(0,1) = w1p-w3+w4p; 00892 myA(0,2) = w1p-w3-w4p; 00893 myA(0,3) = w1-w5; 00894 myA(1,0) = w1-w3p+w4; 00895 myA(1,1) = w1p; 00896 myA(1,2) = w1p-w5p; 00897 myA(1,3) = w1-w3p-w4; 00898 myA(2,0) = w1+w3p+w4; 00899 myA(2,1) = w1p+w5p; 00900 myA(2,2) = w1p; 00901 myA(2,3) = w1+w3p-w4; 00902 myA(3,0) = w1+w5; 00903 myA(3,1) = w1p+w3+w4p; 00904 myA(3,2) = w1p+w3-w4p; 00905 myA(3,3) = w1; 00906 myb(0) = 2*w1; 00907 myb(1) = 2*w1p; 00908 myb(2) = 2*w1p; 00909 myb(3) = 2*w1; 00910 myc(0) = onehalf - w2; 00911 myc(1) = onehalf - w2p; 00912 myc(2) = onehalf + w2p; 00913 myc(3) = onehalf + w2; 00914 this->setMyDescription(myDescription.str()); 00915 this->setMy_A(myA); 00916 this->setMy_b(myb); 00917 this->setMy_c(myc); 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 myDescription; 00932 myDescription << 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 myNumStages = 2; 00943 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00944 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00945 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00946 Scalar one = ST::one(); 00947 Scalar onequarter = as<Scalar>( one/(4*one) ); 00948 Scalar onehalf = as<Scalar>( one/(2*one) ); 00949 myA(0,0) = onequarter; 00950 myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) ); 00951 myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) ); 00952 myA(1,1) = onequarter; 00953 myb(0) = onehalf; 00954 myb(1) = onehalf; 00955 myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) ); 00956 myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) ); 00957 this->setMyDescription(myDescription.str()); 00958 this->setMy_A(myA); 00959 this->setMy_b(myb); 00960 this->setMy_c(myc); 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 myDescription; 00975 myDescription << 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 myNumStages = 1; 00989 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 00990 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 00991 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 00992 Scalar onehalf = ST::one()/(2*ST::one()); 00993 Scalar one = ST::one(); 00994 myA(0,0) = onehalf; 00995 myb(0) = one; 00996 myc(0) = onehalf; 00997 this->setMyDescription(myDescription.str()); 00998 this->setMy_A(myA); 00999 this->setMy_b(myb); 01000 this->setMy_c(myc); 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 myDescription; 01015 myDescription << 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 myNumStages = 2; 01026 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01027 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01028 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = onefourth; 01037 myA(0,1) = onefourth-alpha; 01038 myA(1,0) = onefourth+alpha; 01039 myA(1,1) = onefourth; 01040 myb(0) = onehalf; 01041 myb(1) = onehalf; 01042 myc(0) = onehalf-alpha; 01043 myc(1) = onehalf+alpha; 01044 this->setMyDescription(myDescription.str()); 01045 this->setMy_A(myA); 01046 this->setMy_b(myb); 01047 this->setMy_c(myc); 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 myDescription; 01062 myDescription << 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 myNumStages = 3; 01074 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01075 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01076 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = as<Scalar>(5*one/(36*one)); 01088 myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15; 01089 myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30; 01090 myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24; 01091 myA(1,1) = as<Scalar>(2*one/(9*one)); 01092 myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24; 01093 myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30; 01094 myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15; 01095 myA(2,2) = as<Scalar>(5*one/(36*one)); 01096 myb(0) = as<Scalar>(5*one/(18*one)); 01097 myb(1) = as<Scalar>(4*one/(9*one)); 01098 myb(2) = as<Scalar>(5*one/(18*one)); 01099 myc(0) = as<Scalar>(one/(2*one))-sqrt15over10; 01100 myc(1) = as<Scalar>(one/(2*one)); 01101 myc(2) = as<Scalar>(one/(2*one))+sqrt15over10; 01102 this->setMyDescription(myDescription.str()); 01103 this->setMy_A(myA); 01104 this->setMy_b(myb); 01105 this->setMy_c(myc); 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 myDescription; 01120 myDescription << 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 myNumStages = 1; 01130 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01131 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01132 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01133 Scalar one = ST::one(); 01134 Scalar zero = ST::zero(); 01135 myA(0,0) = one; 01136 myb(0) = one; 01137 myc(0) = zero; 01138 this->setMyDescription(myDescription.str()); 01139 this->setMy_A(myA); 01140 this->setMy_b(myb); 01141 this->setMy_c(myc); 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 myDescription; 01156 myDescription << 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 myNumStages = 2; 01167 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01168 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01169 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01170 Scalar zero = ST::zero(); 01171 Scalar one = ST::one(); 01172 myA(0,0) = as<Scalar>(one/(4*one)); 01173 myA(0,1) = as<Scalar>(-one/(4*one)); 01174 myA(1,0) = as<Scalar>(one/(4*one)); 01175 myA(1,1) = as<Scalar>(5*one/(12*one)); 01176 myb(0) = as<Scalar>(one/(4*one)); 01177 myb(1) = as<Scalar>(3*one/(4*one)); 01178 myc(0) = zero; 01179 myc(1) = as<Scalar>(2*one/(3*one)); 01180 this->setMyDescription(myDescription.str()); 01181 this->setMy_A(myA); 01182 this->setMy_b(myb); 01183 this->setMy_c(myc); 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 myDescription; 01198 myDescription << 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 myNumStages = 3; 01210 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01211 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01212 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01213 Scalar zero = ST::zero(); 01214 Scalar one = ST::one(); 01215 myA(0,0) = as<Scalar>(one/(9*one)); 01216 myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) ); 01217 myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) ); 01218 myA(1,0) = as<Scalar>(one/(9*one)); 01219 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) ); 01220 myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) ); 01221 myA(2,0) = as<Scalar>(one/(9*one)); 01222 myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) ); 01223 myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) ); 01224 myb(0) = as<Scalar>(one/(9*one)); 01225 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) ); 01226 myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) ); 01227 myc(0) = zero; 01228 myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) ); 01229 myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) ); 01230 this->setMyDescription(myDescription.str()); 01231 this->setMy_A(myA); 01232 this->setMy_b(myb); 01233 this->setMy_c(myc); 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 myDescription; 01248 myDescription << 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 myNumStages = 1; 01258 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01259 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01260 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01261 Scalar one = ST::one(); 01262 myA(0,0) = one; 01263 myb(0) = one; 01264 myc(0) = one; 01265 this->setMyDescription(myDescription.str()); 01266 this->setMy_A(myA); 01267 this->setMy_b(myb); 01268 this->setMy_c(myc); 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 myDescription; 01283 myDescription << 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 myNumStages = 2; 01294 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01295 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01296 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01297 Scalar one = ST::one(); 01298 myA(0,0) = as<Scalar>( 5*one/(12*one) ); 01299 myA(0,1) = as<Scalar>( -one/(12*one) ); 01300 myA(1,0) = as<Scalar>( 3*one/(4*one) ); 01301 myA(1,1) = as<Scalar>( one/(4*one) ); 01302 myb(0) = as<Scalar>( 3*one/(4*one) ); 01303 myb(1) = as<Scalar>( one/(4*one) ); 01304 myc(0) = as<Scalar>( one/(3*one) ); 01305 myc(1) = one; 01306 this->setMyDescription(myDescription.str()); 01307 this->setMy_A(myA); 01308 this->setMy_b(myb); 01309 this->setMy_c(myc); 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 myDescription; 01324 myDescription << 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 myNumStages = 3; 01336 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01337 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01338 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01339 Scalar one = ST::one(); 01340 myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) ); 01341 myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) ); 01342 myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) ); 01343 myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) ); 01344 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) ); 01345 myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) ); 01346 myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) ); 01347 myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) ); 01348 myA(2,2) = as<Scalar>( one/(9*one) ); 01349 myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) ); 01350 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) ); 01351 myb(2) = as<Scalar>( one/(9*one) ); 01352 myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) ); 01353 myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) ); 01354 myc(2) = one; 01355 this->setMyDescription(myDescription.str()); 01356 this->setMy_A(myA); 01357 this->setMy_b(myb); 01358 this->setMy_c(myc); 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 myDescription; 01373 myDescription << 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 myNumStages = 2; 01384 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01385 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01386 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01387 Scalar zero = ST::zero(); 01388 Scalar one = ST::one(); 01389 myA(0,0) = zero; 01390 myA(0,1) = zero; 01391 myA(1,0) = as<Scalar>( one/(2*one) ); 01392 myA(1,1) = as<Scalar>( one/(2*one) ); 01393 myb(0) = as<Scalar>( one/(2*one) ); 01394 myb(1) = as<Scalar>( one/(2*one) ); 01395 myc(0) = zero; 01396 myc(1) = one; 01397 this->setMyDescription(myDescription.str()); 01398 this->setMy_A(myA); 01399 this->setMy_b(myb); 01400 this->setMy_c(myc); 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 myDescription; 01415 myDescription << 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 myNumStages = 3; 01427 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01428 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01429 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01430 Scalar zero = ST::zero(); 01431 Scalar one = ST::one(); 01432 myA(0,0) = zero; 01433 myA(0,1) = zero; 01434 myA(0,2) = zero; 01435 myA(1,0) = as<Scalar>( (5*one)/(24*one) ); 01436 myA(1,1) = as<Scalar>( (one)/(3*one) ); 01437 myA(1,2) = as<Scalar>( (-one)/(24*one) ); 01438 myA(2,0) = as<Scalar>( (one)/(6*one) ); 01439 myA(2,1) = as<Scalar>( (2*one)/(3*one) ); 01440 myA(2,2) = as<Scalar>( (1*one)/(6*one) ); 01441 myb(0) = as<Scalar>( (one)/(6*one) ); 01442 myb(1) = as<Scalar>( (2*one)/(3*one) ); 01443 myb(2) = as<Scalar>( (1*one)/(6*one) ); 01444 myc(0) = zero; 01445 myc(1) = as<Scalar>( one/(2*one) ); 01446 myc(2) = one; 01447 this->setMyDescription(myDescription.str()); 01448 this->setMy_A(myA); 01449 this->setMy_b(myb); 01450 this->setMy_c(myc); 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 myDescription; 01468 myDescription << 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 myNumStages = 4; 01481 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01482 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01483 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01484 Scalar zero = ST::zero(); 01485 Scalar one = ST::one(); 01486 myA(0,0) = zero; 01487 myA(0,1) = zero; 01488 myA(0,2) = zero; 01489 myA(0,3) = zero; 01490 myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) ); 01491 myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) ); 01492 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) ); 01493 myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) ); 01494 myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) ); 01495 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) ); 01496 myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) ); 01497 myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) ); 01498 myA(3,0) = as<Scalar>( one/(12*one) ); 01499 myA(3,1) = as<Scalar>( 5*one/(12*one) ); 01500 myA(3,2) = as<Scalar>( 5*one/(12*one) ); 01501 myA(3,3) = as<Scalar>( one/(12*one) ); 01502 myb(0) = as<Scalar>( one/(12*one) ); 01503 myb(1) = as<Scalar>( 5*one/(12*one) ); 01504 myb(2) = as<Scalar>( 5*one/(12*one) ); 01505 myb(3) = as<Scalar>( one/(12*one) ); 01506 myc(0) = zero; 01507 myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) ); 01508 myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) ); 01509 myc(3) = one; 01510 this->setMyDescription(myDescription.str()); 01511 this->setMy_A(myA); 01512 this->setMy_b(myb); 01513 this->setMy_c(myc); 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 myDescription; 01528 myDescription << 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 myNumStages = 2; 01539 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01540 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01541 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01542 Scalar zero = ST::zero(); 01543 Scalar one = ST::one(); 01544 myA(0,0) = as<Scalar>( one/(2*one) ); 01545 myA(0,1) = zero; 01546 myA(1,0) = as<Scalar>( one/(2*one) ); 01547 myA(1,1) = zero; 01548 myb(0) = as<Scalar>( one/(2*one) ); 01549 myb(1) = as<Scalar>( one/(2*one) ); 01550 myc(0) = zero; 01551 myc(1) = one; 01552 this->setMyDescription(myDescription.str()); 01553 this->setMy_A(myA); 01554 this->setMy_b(myb); 01555 this->setMy_c(myc); 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 myDescription; 01570 myDescription << 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 myNumStages = 3; 01582 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01583 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01584 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01585 Scalar zero = ST::zero(); 01586 Scalar one = ST::one(); 01587 myA(0,0) = as<Scalar>( one/(6*one) ); 01588 myA(0,1) = as<Scalar>( -one/(6*one) ); 01589 myA(0,2) = zero; 01590 myA(1,0) = as<Scalar>( one/(6*one) ); 01591 myA(1,1) = as<Scalar>( one/(3*one) ); 01592 myA(1,2) = zero; 01593 myA(2,0) = as<Scalar>( one/(6*one) ); 01594 myA(2,1) = as<Scalar>( 5*one/(6*one) ); 01595 myA(2,2) = zero; 01596 myb(0) = as<Scalar>( one/(6*one) ); 01597 myb(1) = as<Scalar>( 2*one/(3*one) ); 01598 myb(2) = as<Scalar>( one/(6*one) ); 01599 myc(0) = zero; 01600 myc(1) = as<Scalar>( one/(2*one) ); 01601 myc(2) = one; 01602 this->setMyDescription(myDescription.str()); 01603 this->setMy_A(myA); 01604 this->setMy_b(myb); 01605 this->setMy_c(myc); 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 myDescription; 01620 myDescription << 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 myNumStages = 4; 01633 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01634 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01635 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01636 Scalar zero = ST::zero(); 01637 Scalar one = ST::one(); 01638 myA(0,0) = as<Scalar>( one/(12*one) ); 01639 myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) ); 01640 myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) ); 01641 myA(0,3) = zero; 01642 myA(1,0) = as<Scalar>( one/(12*one) ); 01643 myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) ); 01644 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) ); 01645 myA(1,3) = zero; 01646 myA(2,0) = as<Scalar>( one/(12*one) ); 01647 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) ); 01648 myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) ); 01649 myA(2,3) = zero; 01650 myA(3,0) = as<Scalar>( one/(12*one) ); 01651 myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) ); 01652 myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) ); 01653 myA(3,3) = zero; 01654 myb(0) = as<Scalar>( one/(12*one) ); 01655 myb(1) = as<Scalar>( 5*one/(12*one) ); 01656 myb(2) = as<Scalar>( 5*one/(12*one) ); 01657 myb(3) = as<Scalar>( one/(12*one) ); 01658 myc(0) = zero; 01659 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) ); 01660 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) ); 01661 myc(3) = one; 01662 this->setMyDescription(myDescription.str()); 01663 this->setMy_A(myA); 01664 this->setMy_b(myb); 01665 this->setMy_c(myc); 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 myDescription; 01680 myDescription << 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 myNumStages = 2; 01691 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01692 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01693 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01694 Scalar zero = ST::zero(); 01695 Scalar one = ST::one(); 01696 myA(0,0) = as<Scalar>( one/(2*one) ); 01697 myA(0,1) = as<Scalar>( -one/(2*one) ); 01698 myA(1,0) = as<Scalar>( one/(2*one) ); 01699 myA(1,1) = as<Scalar>( one/(2*one) ); 01700 myb(0) = as<Scalar>( one/(2*one) ); 01701 myb(1) = as<Scalar>( one/(2*one) ); 01702 myc(0) = zero; 01703 myc(1) = one; 01704 this->setMyDescription(myDescription.str()); 01705 this->setMy_A(myA); 01706 this->setMy_b(myb); 01707 this->setMy_c(myc); 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 myDescription; 01722 myDescription << 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 myNumStages = 3; 01734 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01735 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01736 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01737 Scalar zero = ST::zero(); 01738 Scalar one = ST::one(); 01739 myA(0,0) = as<Scalar>( one/(6*one) ); 01740 myA(0,1) = as<Scalar>( -one/(3*one) ); 01741 myA(0,2) = as<Scalar>( one/(6*one) ); 01742 myA(1,0) = as<Scalar>( one/(6*one) ); 01743 myA(1,1) = as<Scalar>( 5*one/(12*one) ); 01744 myA(1,2) = as<Scalar>( -one/(12*one) ); 01745 myA(2,0) = as<Scalar>( one/(6*one) ); 01746 myA(2,1) = as<Scalar>( 2*one/(3*one) ); 01747 myA(2,2) = as<Scalar>( one/(6*one) ); 01748 myb(0) = as<Scalar>( one/(6*one) ); 01749 myb(1) = as<Scalar>( 2*one/(3*one) ); 01750 myb(2) = as<Scalar>( one/(6*one) ); 01751 myc(0) = zero; 01752 myc(1) = as<Scalar>( one/(2*one) ); 01753 myc(2) = one; 01754 this->setMyDescription(myDescription.str()); 01755 this->setMy_A(myA); 01756 this->setMy_b(myb); 01757 this->setMy_c(myc); 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 myDescription; 01772 myDescription << 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 myNumStages = 4; 01785 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01786 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01787 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01788 Scalar zero = ST::zero(); 01789 Scalar one = ST::one(); 01790 myA(0,0) = as<Scalar>( one/(12*one) ); 01791 myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) ); 01792 myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) ); 01793 myA(0,3) = as<Scalar>( -one/(12*one) ); 01794 myA(1,0) = as<Scalar>( one/(12*one) ); 01795 myA(1,1) = as<Scalar>( one/(4*one) ); 01796 myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) ); 01797 myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) ); 01798 myA(2,0) = as<Scalar>( one/(12*one) ); 01799 myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) ); 01800 myA(2,2) = as<Scalar>( one/(4*one) ); 01801 myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) ); 01802 myA(3,0) = as<Scalar>( one/(12*one) ); 01803 myA(3,1) = as<Scalar>( 5*one/(12*one) ); 01804 myA(3,2) = as<Scalar>( 5*one/(12*one) ); 01805 myA(3,3) = as<Scalar>( one/(12*one) ); 01806 myb(0) = as<Scalar>( one/(12*one) ); 01807 myb(1) = as<Scalar>( 5*one/(12*one) ); 01808 myb(2) = as<Scalar>( 5*one/(12*one) ); 01809 myb(3) = as<Scalar>( one/(12*one) ); 01810 myc(0) = zero; 01811 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) ); 01812 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) ); 01813 myc(3) = one; 01814 this->setMyDescription(myDescription.str()); 01815 this->setMy_A(myA); 01816 this->setMy_b(myb); 01817 this->setMy_c(myc); 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 myDescription; 01833 myDescription << 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 myNumStages = 5; 01847 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01848 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01849 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 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 myA(0,0) = gamma; 01855 myA(0,1) = zero; 01856 myA(0,2) = zero; 01857 myA(0,3) = zero; 01858 myA(0,4) = zero; 01859 01860 myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) ); 01861 myA(1,1) = gamma; 01862 myA(1,2) = zero; 01863 myA(1,3) = zero; 01864 myA(1,4) = zero; 01865 01866 myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) ); 01867 myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) ); 01868 myA(2,2) = gamma; 01869 myA(2,3) = zero; 01870 myA(2,4) = zero; 01871 01872 myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) ); 01873 myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) ); 01874 myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) ); 01875 myA(3,3) = gamma; 01876 myA(3,4) = zero; 01877 01878 myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) ); 01879 myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) ); 01880 myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) ); 01881 myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) ); 01882 myA(4,4) = gamma; 01883 01884 myb(0) = zero; 01885 myb(1) = zero; 01886 myb(2) = as<Scalar>( one/(9*one) ); 01887 myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) ); 01888 myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) ); 01889 01890 myc(0) = gamma; 01891 myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) ); 01892 myc(2) = one; 01893 myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) ); 01894 myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) ); 01895 01896 this->setMyDescription(myDescription.str()); 01897 this->setMy_A(myA); 01898 this->setMy_b(myb); 01899 this->setMy_c(myc); 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 myDescription; 01914 myDescription << 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 myNumStages = 5; 01929 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 01930 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 01931 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 01932 Scalar zero = ST::zero(); 01933 Scalar one = ST::one(); 01934 Scalar onequarter = as<Scalar>( one/(4*one) ); 01935 myA(0,0) = onequarter; 01936 myA(0,1) = zero; 01937 myA(0,2) = zero; 01938 myA(0,3) = zero; 01939 myA(0,4) = zero; 01940 01941 myA(1,0) = as<Scalar>( one / (2*one) ); 01942 myA(1,1) = onequarter; 01943 myA(1,2) = zero; 01944 myA(1,3) = zero; 01945 myA(1,4) = zero; 01946 01947 myA(2,0) = as<Scalar>( 17*one/(50*one) ); 01948 myA(2,1) = as<Scalar>( -one/(25*one) ); 01949 myA(2,2) = onequarter; 01950 myA(2,3) = zero; 01951 myA(2,4) = zero; 01952 01953 myA(3,0) = as<Scalar>( 371*one/(1360*one) ); 01954 myA(3,1) = as<Scalar>( -137*one/(2720*one) ); 01955 myA(3,2) = as<Scalar>( 15*one/(544*one) ); 01956 myA(3,3) = onequarter; 01957 myA(3,4) = zero; 01958 01959 myA(4,0) = as<Scalar>( 25*one/(24*one) ); 01960 myA(4,1) = as<Scalar>( -49*one/(48*one) ); 01961 myA(4,2) = as<Scalar>( 125*one/(16*one) ); 01962 myA(4,3) = as<Scalar>( -85*one/(12*one) ); 01963 myA(4,4) = onequarter; 01964 01965 myb(0) = as<Scalar>( 25*one/(24*one) ); 01966 myb(1) = as<Scalar>( -49*one/(48*one) ); 01967 myb(2) = as<Scalar>( 125*one/(16*one) ); 01968 myb(3) = as<Scalar>( -85*one/(12*one) ); 01969 myb(4) = onequarter; 01970 01971 /* 01972 // Alternate version 01973 myb(0) = as<Scalar>( 59*one/(48*one) ); 01974 myb(1) = as<Scalar>( -17*one/(96*one) ); 01975 myb(2) = as<Scalar>( 225*one/(32*one) ); 01976 myb(3) = as<Scalar>( -85*one/(12*one) ); 01977 myb(4) = zero; 01978 */ 01979 myc(0) = onequarter; 01980 myc(1) = as<Scalar>( 3*one/(4*one) ); 01981 myc(2) = as<Scalar>( 11*one/(20*one) ); 01982 myc(3) = as<Scalar>( one/(2*one) ); 01983 myc(4) = one; 01984 01985 this->setMyDescription(myDescription.str()); 01986 this->setMy_A(myA); 01987 this->setMy_b(myb); 01988 this->setMy_c(myc); 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 myDescription; 02003 myDescription << 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 myNumStages = 3; 02017 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages); 02018 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages); 02019 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages); 02020 Scalar zero = ST::zero(); 02021 Scalar one = ST::one(); 02022 Scalar pi = as<Scalar>(4*one)*std::atan(one); 02023 Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) ); 02024 Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) ); 02025 myA(0,0) = gamma; 02026 myA(0,1) = zero; 02027 myA(0,2) = zero; 02028 02029 myA(1,0) = as<Scalar>( one/(2*one) - gamma ); 02030 myA(1,1) = gamma; 02031 myA(1,2) = zero; 02032 02033 myA(2,0) = as<Scalar>( 2*gamma ); 02034 myA(2,1) = as<Scalar>( one - 4*gamma ); 02035 myA(2,2) = gamma; 02036 02037 myb(0) = delta; 02038 myb(1) = as<Scalar>( one-2*delta ); 02039 myb(2) = delta; 02040 02041 myc(0) = gamma; 02042 myc(1) = as<Scalar>( one/(2*one) ); 02043 myc(2) = as<Scalar>( one - gamma ); 02044 02045 this->setMyDescription(myDescription.str()); 02046 this->setMy_A(myA); 02047 this->setMy_b(myb); 02048 this->setMy_c(myc); 02049 this->setMy_order(4); 02050 } 02051 }; 02052 02053 02054 } // namespace Rythmos 02055 02056 02057 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP
1.7.4
