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