1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2011 Fabien Le Floc'h 5 6 This file is part of QuantLib, a free-software/open-source library 7 for financial quantitative analysts and developers - http://quantlib.org/ 8 9 QuantLib is free software: you can redistribute it and/or modify it 10 under the terms of the QuantLib license. You should have received a 11 copy of the license along with this program; if not, please email 12 <quantlib-dev@lists.sf.net>. The license is also available online at 13 <http://quantlib.org/license.shtml>. 14 15 This program is distributed in the hope that it will be useful, but WITHOUT 16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 17 FOR A PARTICULAR PURPOSE. See the license for more details. 18 */ 19 20 /*! \file trbdf2.hpp 21 \brief TR-BDF2 scheme for finite difference methods 22 */ 23 24 #ifndef quantlib_trbdf2_hpp 25 #define quantlib_trbdf2_hpp 26 27 #include <ql/methods/finitedifferences/finitedifferencemodel.hpp> 28 29 namespace QuantLib { 30 31 //! TR-BDF2 scheme for finite difference methods 32 /*! See <http://ssrn.com/abstract=1648878> for details. 33 34 In this implementation, the passed operator must be derived 35 from either TimeConstantOperator or TimeDependentOperator. 36 Also, it must implement at least the following interface: 37 38 \code 39 typedef ... array_type; 40 41 // copy constructor/assignment 42 // (these will be provided by the compiler if none is defined) 43 Operator(const Operator&); 44 Operator& operator=(const Operator&); 45 46 // inspectors 47 Size size(); 48 49 // modifiers 50 void setTime(Time t); 51 52 // operator interface 53 array_type applyTo(const array_type&); 54 array_type solveFor(const array_type&); 55 static Operator identity(Size size); 56 57 // operator algebra 58 Operator operator*(Real, const Operator&); 59 Operator operator+(const Operator&, const Operator&); 60 Operator operator+(const Operator&, const Operator&); 61 \endcode 62 63 \warning The differential operator must be linear for 64 this evolver to work. 65 66 \ingroup findiff 67 */ 68 69 // NOTE: There is room for performance improvement especially in 70 // the array manipulation 71 72 template <class Operator> 73 class TRBDF2 { 74 public: 75 // typedefs 76 typedef OperatorTraits<Operator> traits; 77 typedef typename traits::operator_type operator_type; 78 typedef typename traits::array_type array_type; 79 typedef typename traits::bc_set bc_set; 80 typedef typename traits::condition_type condition_type; 81 // constructors TRBDF2(const operator_type & L,const bc_set & bcs)82 TRBDF2(const operator_type& L, 83 const bc_set& bcs) 84 : L_(L), I_(operator_type::identity(L.size())), 85 dt_(0.0), bcs_(bcs), alpha_(2.0-sqrt(2.0)) {} 86 void step(array_type& a, 87 Time t); setStep(Time dt)88 void setStep(Time dt) { 89 dt_ = dt; 90 91 implicitPart_ = I_ + 0.5*alpha_*dt_*L_; 92 explicitTrapezoidalPart_ = I_ - 0.5*alpha_*dt_*L_; 93 explicitBDF2PartFull_ = 94 -(1.0-alpha_)*(1.0-alpha_)/(alpha_*(2.0-alpha_))*I_; 95 explicitBDF2PartMid_ = 1.0/(alpha_*(2.0-alpha_))*I_; 96 } 97 private: 98 Real alpha_; 99 operator_type L_, I_, explicitTrapezoidalPart_, 100 explicitBDF2PartFull_,explicitBDF2PartMid_, implicitPart_; 101 Time dt_; 102 bc_set bcs_; 103 array_type aInit_; 104 }; 105 106 107 // inline definitions 108 109 template <class Operator> step(array_type & a,Time t)110 inline void TRBDF2<Operator>::step(array_type& a, Time t) { 111 Size i; 112 Array aInit(a.size()); 113 for (i=0; i<a.size();i++) { 114 aInit[i] = a[i]; 115 } 116 aInit_ = aInit; 117 for (i=0; i<bcs_.size(); i++) 118 bcs_[i]->setTime(t); 119 //trapezoidal explicit part 120 if (L_.isTimeDependent()) { 121 L_.setTime(t); 122 explicitTrapezoidalPart_ = I_ - 0.5*alpha_*dt_*L_; 123 } 124 for (i=0; i<bcs_.size(); i++) 125 bcs_[i]->applyBeforeApplying(explicitTrapezoidalPart_); 126 a = explicitTrapezoidalPart_.applyTo(a); 127 for (i=0; i<bcs_.size(); i++) 128 bcs_[i]->applyAfterApplying(a); 129 130 // trapezoidal implicit part 131 if (L_.isTimeDependent()) { 132 L_.setTime(t-dt_); 133 implicitPart_ = I_ + 0.5*alpha_*dt_*L_; 134 } 135 for (i=0; i<bcs_.size(); i++) 136 bcs_[i]->applyBeforeSolving(implicitPart_,a); 137 a = implicitPart_.solveFor(a); 138 for (i=0; i<bcs_.size(); i++) 139 bcs_[i]->applyAfterSolving(a); 140 141 142 // BDF2 explicit part 143 if (L_.isTimeDependent()) { 144 L_.setTime(t); 145 } 146 for (i=0; i<bcs_.size(); i++) { 147 bcs_[i]->applyBeforeApplying(explicitBDF2PartFull_); 148 } 149 array_type b0 = explicitBDF2PartFull_.applyTo(aInit_); 150 for (i=0; i<bcs_.size(); i++) 151 bcs_[i]->applyAfterApplying(b0); 152 153 for (i=0; i<bcs_.size(); i++) { 154 bcs_[i]->applyBeforeApplying(explicitBDF2PartMid_); 155 } 156 array_type b1 = explicitBDF2PartMid_.applyTo(a); 157 for (i=0; i<bcs_.size(); i++) 158 bcs_[i]->applyAfterApplying(b1); 159 a = b0+b1; 160 161 // reuse implicit part - works only for alpha=2-sqrt(2) 162 for (i=0; i<bcs_.size(); i++) 163 bcs_[i]->applyBeforeSolving(implicitPart_,a); 164 a = implicitPart_.solveFor(a); 165 for (i=0; i<bcs_.size(); i++) 166 bcs_[i]->applyAfterSolving(a); 167 168 } 169 170 } 171 172 #endif 173