QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
trbdf2.hpp
Go to the documentation of this file.
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
28
29namespace 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
79 typedef typename traits::bc_set bc_set;
81 // constructors
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);
88 void setStep(Time dt) {
89 dt_ = dt;
90
94 -(1.0-alpha_)*(1.0-alpha_)/(alpha_*(2.0-alpha_))*I_;
96 }
97 private:
104 };
105
106
107 // inline definitions
108
109 template <class Operator>
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
1-D array used in linear algebra.
Definition: array.hpp:52
std::vector< ext::shared_ptr< bc_type > > bc_set
Operator::array_type array_type
condition to be applied at every time step
TR-BDF2 scheme for finite difference methods.
Definition: trbdf2.hpp:73
traits::operator_type operator_type
Definition: trbdf2.hpp:77
operator_type explicitBDF2PartFull_
Definition: trbdf2.hpp:100
array_type aInit_
Definition: trbdf2.hpp:103
operator_type explicitBDF2PartMid_
Definition: trbdf2.hpp:100
operator_type explicitTrapezoidalPart_
Definition: trbdf2.hpp:99
TRBDF2(const operator_type &L, const bc_set &bcs)
Definition: trbdf2.hpp:82
operator_type implicitPart_
Definition: trbdf2.hpp:100
void setStep(Time dt)
Definition: trbdf2.hpp:88
operator_type L_
Definition: trbdf2.hpp:99
traits::bc_set bc_set
Definition: trbdf2.hpp:79
operator_type I_
Definition: trbdf2.hpp:99
traits::condition_type condition_type
Definition: trbdf2.hpp:80
OperatorTraits< Operator > traits
Definition: trbdf2.hpp:76
traits::array_type array_type
Definition: trbdf2.hpp:78
void step(array_type &a, Time t)
Definition: trbdf2.hpp:110
const DefaultType & t
generic finite difference model
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35