1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
Copyright (C) 2008 Andreas Gaida
Copyright (C) 2008 Ralph Schreyer
Copyright (C) 2008 Klaus Spanderen
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
<quantlib-dev@lists.sf.net>. The license is also available online at
<http://quantlib.org/license.shtml>.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
#include <ql/methods/finitedifferences/operators/secondordermixedderivativeop.hpp>
namespace QuantLib {
SecondOrderMixedDerivativeOp::SecondOrderMixedDerivativeOp(
Size d0, Size d1,
const ext::shared_ptr<FdmMesher>& mesher)
: NinePointLinearOp(d0, d1, mesher) {
const ext::shared_ptr<FdmLinearOpLayout> layout = mesher->layout();
const FdmLinearOpIterator endIter = layout->end();
for (FdmLinearOpIterator iter = layout->begin(); iter!=endIter; ++iter) {
const Size i = iter.index();
const Real hm_d0 = mesher->dminus(iter, d0_);
const Real hp_d0 = mesher->dplus(iter, d0_);
const Real hm_d1 = mesher->dminus(iter, d1_);
const Real hp_d1 = mesher->dplus(iter, d1_);
const Real zetam1 = hm_d0*(hm_d0+hp_d0);
const Real zeta0 = hm_d0*hp_d0;
const Real zetap1 = hp_d0*(hm_d0+hp_d0);
const Real phim1 = hm_d1*(hm_d1+hp_d1);
const Real phi0 = hm_d1*hp_d1;
const Real phip1 = hp_d1*(hm_d1+hp_d1);
const Size c0 = iter.coordinates()[d0_];
const Size c1 = iter.coordinates()[d1_];
if (c0 == 0 && c1 == 0) {
// lower left corner
a00_[i] = a01_[i] = a02_[i] = a10_[i] = a20_[i] = 0.0;
a21_[i] = a12_[i] = -(a11_[i] = a22_[i] = 1.0/(hp_d0*hp_d1));
}
else if (c0 == layout->dim()[d0_]-1 && c1 == 0) {
// upper left corner
a22_[i] = a21_[i] = a20_[i] = a10_[i] = a00_[i] = 0.0;
a11_[i] = a02_[i] = -(a01_[i] = a12_[i] = 1.0/(hm_d0*hp_d1));
}
else if (c0 == 0 && c1 == layout->dim()[d1_]-1) {
// lower right corner
a00_[i] = a01_[i] = a02_[i] = a12_[i] = a22_[i] = 0.0;
a20_[i] = a11_[i] = -(a10_[i] = a21_[i] = 1.0/(hp_d0*hm_d1));
}
else if (c0 == layout->dim()[d0_]-1 && c1 == layout->dim()[d1_]-1) {
// upper right corner
a20_[i] = a21_[i] = a22_[i] = a12_[i] = a02_[i] = 0.0;
a10_[i] = a01_[i] = -(a00_[i] = a11_[i] = 1.0/(hm_d0*hm_d1));
}
else if (c0 == 0) {
// lower side
a00_[i] = a01_[i] = a02_[i] = 0.0;
a20_[i] = -(a10_[i] = hp_d1/(hp_d0*phim1));
a11_[i] = -(a21_[i] = (hp_d1-hm_d1)/(hp_d0*phi0));
a12_[i] = -(a22_[i] = hm_d1/(hp_d0*phip1));
}
else if (c0 == layout->dim()[d0_]-1) {
// upper side
a20_[i] = a21_[i] = a22_[i] = 0.0;
a10_[i] = -(a00_[i] = hp_d1/(hm_d0*phim1));
a01_[i] = -(a11_[i] = (hp_d1-hm_d1)/(hm_d0*phi0));
a02_[i] = -(a12_[i] = hm_d1/(hm_d0*phip1));
}
else if (c1 == 0) {
// left side
a00_[i] = a10_[i] = a20_[i] = 0.0;
a02_[i] = -(a01_[i] = hp_d0/(zetam1*hp_d1));
a11_[i] = -(a12_[i] = (hp_d0-hm_d0)/(zeta0*hp_d1));
a21_[i] = -(a22_[i] = hm_d0/(zetap1*hp_d1));
}
else if (c1 == layout->dim()[d1_]-1) {
// right side
a22_[i] = a12_[i] = a02_[i] = 0.0;
a01_[i] = -(a00_[i] = hp_d0/(zetam1*hm_d1));
a10_[i] = -(a11_[i] = (hp_d0-hm_d0)/(zeta0*hm_d1));
a20_[i] = -(a21_[i] = hm_d0/(zetap1*hm_d1));
}
else {
a00_[i] = hp_d0*hp_d1/(zetam1*phim1);
a10_[i] = -(hp_d0-hm_d0)*hp_d1/(zeta0*phim1);
a20_[i] = -hm_d0*hp_d1/(zetap1*phim1);
a01_[i] = -hp_d0*(hp_d1-hm_d1)/(zetam1*phi0);
a11_[i] = (hp_d0-hm_d0)*(hp_d1-hm_d1)/(zeta0*phi0);
a21_[i] = hm_d0*(hp_d1-hm_d1)/(zetap1*phi0);
a02_[i] = -hp_d0*hm_d1/(zetam1*phip1);
a12_[i] = hm_d1*(hp_d0-hm_d0)/(zeta0*phip1);
a22_[i] = hm_d0*hm_d1/(zetap1*phip1);
}
}
}
}
|
