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