1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2002, 2003 Ferdinando Ametrano 5 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl 6 7 This file is part of QuantLib, a free-software/open-source library 8 for financial quantitative analysts and developers - http://quantlib.org/ 9 10 QuantLib is free software: you can redistribute it and/or modify it 11 under the terms of the QuantLib license. You should have received a 12 copy of the license along with this program; if not, please email 13 <quantlib-dev@lists.sf.net>. The license is also available online at 14 <http://quantlib.org/license.shtml>. 15 16 This program is distributed in the hope that it will be useful, but WITHOUT 17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 18 FOR A PARTICULAR PURPOSE. See the license for more details. 19 */ 20 21 /*! \file mixedscheme.hpp 22 \brief Mixed (explicit/implicit) scheme for finite difference methods 23 */ 24 25 #ifndef quantlib_mixed_scheme_hpp 26 #define quantlib_mixed_scheme_hpp 27 28 #include <ql/methods/finitedifferences/finitedifferencemodel.hpp> 29 30 namespace QuantLib { 31 32 //! Mixed (explicit/implicit) scheme for finite difference methods 33 /*! In this implementation, the passed operator must be derived 34 from either TimeConstantOperator or TimeDependentOperator. 35 Also, it must implement at least the following interface: 36 37 \code 38 typedef ... array_type; 39 40 // copy constructor/assignment 41 // (these will be provided by the compiler if none is defined) 42 Operator(const Operator&); 43 Operator& operator=(const Operator&); 44 45 // inspectors 46 Size size(); 47 48 // modifiers 49 void setTime(Time t); 50 51 // operator interface 52 array_type applyTo(const array_type&); 53 array_type solveFor(const array_type&); 54 static Operator identity(Size size); 55 56 // operator algebra 57 Operator operator*(Real, const Operator&); 58 Operator operator+(const Operator&, const Operator&); 59 Operator operator+(const Operator&, const Operator&); 60 \endcode 61 62 \warning The differential operator must be linear for 63 this evolver to work. 64 65 \todo 66 - derive variable theta schemes 67 - introduce multi time-level schemes. 68 69 \ingroup findiff 70 */ 71 template <class Operator> 72 class MixedScheme { 73 public: 74 // typedefs 75 typedef OperatorTraits<Operator> traits; 76 typedef typename traits::operator_type operator_type; 77 typedef typename traits::array_type array_type; 78 typedef typename traits::bc_set bc_set; 79 typedef typename traits::condition_type condition_type; 80 // constructors MixedScheme(const operator_type & L,Real theta,const bc_set & bcs)81 MixedScheme(const operator_type& L, 82 Real theta, 83 const bc_set& bcs) 84 : L_(L), I_(operator_type::identity(L.size())), 85 dt_(0.0), theta_(theta) , bcs_(bcs) {} 86 void step(array_type& a, 87 Time t); setStep(Time dt)88 void setStep(Time dt) { 89 dt_ = dt; 90 if (theta_!=1.0) // there is an explicit part 91 explicitPart_ = I_-((1.0-theta_) * dt_)*L_; 92 if (theta_!=0.0) // there is an implicit part 93 implicitPart_ = I_+(theta_ * dt_)*L_; 94 } 95 protected: 96 operator_type L_, I_, explicitPart_, implicitPart_; 97 Time dt_; 98 Real theta_; 99 bc_set bcs_; 100 }; 101 102 103 // inline definitions 104 105 template <class Operator> step(array_type & a,Time t)106 inline void MixedScheme<Operator>::step(array_type& a, Time t) { 107 Size i; 108 for (i=0; i<bcs_.size(); i++) 109 bcs_[i]->setTime(t); 110 if (theta_!=1.0) { // there is an explicit part 111 if (L_.isTimeDependent()) { 112 L_.setTime(t); 113 explicitPart_ = I_-((1.0-theta_) * dt_)*L_; 114 } 115 for (i=0; i<bcs_.size(); i++) 116 bcs_[i]->applyBeforeApplying(explicitPart_); 117 a = explicitPart_.applyTo(a); 118 for (i=0; i<bcs_.size(); i++) 119 bcs_[i]->applyAfterApplying(a); 120 } 121 if (theta_!=0.0) { // there is an implicit part 122 if (L_.isTimeDependent()) { 123 L_.setTime(t-dt_); 124 implicitPart_ = I_+(theta_ * dt_)*L_; 125 } 126 for (i=0; i<bcs_.size(); i++) 127 bcs_[i]->applyBeforeSolving(implicitPart_,a); 128 implicitPart_.solveFor(a, a); 129 for (i=0; i<bcs_.size(); i++) 130 bcs_[i]->applyAfterSolving(a); 131 } 132 } 133 134 } 135 136 137 #endif 138