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