1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2012 Peter Caspers
5 Copyright (C) 2013 Klaus Spanderen
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#include "ode.hpp"
22#include "utilities.hpp"
23#include <ql/math/matrixutilities/expm.hpp>
24#include <ql/math/ode/adaptiverungekutta.hpp>
25#include <complex>
26
27
28using namespace QuantLib;
29using namespace boost::unit_test_framework;
30
31using std::exp;
32using std::sin;
33
34namespace {
35
36 struct ode1 {
37 Real operator()(Real x, Real y) const { return y; }
38 };
39
40 struct ode2 {
41 std::complex<Real> operator()(Real x,
42 const std::complex<Real>& y) {
43 return std::complex<Real>(0.0,1.0)*y;
44 }
45 };
46
47 struct ode3 {
48 std::vector<Real> operator()(Real x, const std::vector<Real>& y) {
49 std::vector<Real> r(2);
50 r[0] = y[1]; r[1] = -y[0];
51 return r;
52 }
53 };
54
55 struct ode4 {
56 std::vector<std::complex<Real> > operator()(
57 const std::complex<Real>& x,
58 const std::vector<std::complex<Real> >& y) {
59 std::vector<std::complex<Real> > r(2);
60 r[0] = y[1]; r[1] = -y[0];
61 return r;
62 }
63 };
64
65}
66
67void OdeTest::testAdaptiveRungeKutta() {
68
69 BOOST_TEST_MESSAGE("Testing adaptive Runge Kutta...");
70
71 AdaptiveRungeKutta<Real> rk_real(1E-12,1E-4,0.0);
72 AdaptiveRungeKutta<std::complex<Real> > rk_complex(1E-12,1E-4,0.0);
73 Real tol1 = 5E-10, tol2 = 2E-12, tol3 = 2E-12, tol4 = 2E-12;
74
75 // f'=f, f(0)=1
76 AdaptiveRungeKutta<Real>::OdeFct1d ode1_ = ode1();
77 Real y10=1;
78
79 // f'=f, f(0)=i
80 AdaptiveRungeKutta<std::complex<Real> >::OdeFct1d ode2_ = ode2();
81 std::complex<Real> y20(0.0,1.0);
82
83 // f''=-f, f(0)=0, f'(0)=1
84 AdaptiveRungeKutta<Real>::OdeFct ode3_ = ode3();
85 std::vector<Real> y30(2); y30[0] = 0.0; y30[1] = 1.0;
86
87 // f''=-f, f(0)=1, f'(0)=i
88 AdaptiveRungeKutta<std::complex<Real> >::OdeFct ode4_ = ode4();
89 std::vector<std::complex<Real> > y40(2);
90 y40[0] = 1.0;
91 y40[1] = std::complex<Real>(0.0,1.0);
92
93 Real x=0.0;
94 Real y1 = y10;
95 std::complex<Real> y2 = y20;
96 std::vector<Real> y3 = y30;
97 std::vector<std::complex<Real> > y4 = y40;
98
99 while (x<5.0) {
100 Real exact1 = exp(x: x);
101 std::complex<Real> exact2 =
102 std::exp(z: std::complex<Real>(0.0,x)) * std::complex<Real>(0.0,1.0);
103 Real exact3 = sin(x: x);
104 std::complex<Real> exact4 = std::exp(z: std::complex<Real>(0.0,x));
105
106 if ( std::fabs( x: exact1 - y1 ) > tol1 )
107 BOOST_FAIL("Error in ode #1: exact solution at x=" << x
108 << " is " << exact1
109 << ", numerical solution is " << y1
110 << " difference " << std::fabs(exact1-y1)
111 << " outside tolerance " << tol1);
112 if ( abs( z: exact2 - y2 ) > tol2 )
113 BOOST_FAIL("Error in ode #2: exact solution at x=" << x
114 << " is " << exact2
115 << ", numerical solution is " << y2
116 << " difference " << abs(exact2-y2)
117 << " outside tolerance " << tol2);
118 if ( std::fabs( x: exact3 - y3[0] ) > tol3 )
119 BOOST_FAIL("Error in ode #3: exact solution at x=" << x
120 << " is " << exact3
121 << ", numerical solution is " << y3[0]
122 << " difference " << std::fabs(exact3-y3[0])
123 << " outside tolerance " << tol3);
124 if ( abs( z: exact4 - y4[0] ) > tol4 )
125 BOOST_FAIL("Error in ode #4: exact solution at x=" << x
126 << " is " << exact4
127 << ", numerical solution is " << y4[0]
128 << " difference " << abs(exact4-y4[0])
129 << " outside tolerance " << tol4);
130 x+=0.01;
131 y1=rk_real(ode1_,y10,0.0,x);
132 y2=rk_complex(ode2_,y20,0.0,x);
133 y3=rk_real(ode3_,y30,0.0,x);
134 y4=rk_complex(ode4_,y40,0.0,x);
135 }
136}
137
138namespace {
139 Real frobenuiusNorm(const Matrix& m) {
140 return std::sqrt(x: DotProduct(v1: (m*transpose(m)).diagonal(),
141 v2: Array(m.rows(), 1.0)));
142 }
143}
144
145void OdeTest::testMatrixExponential() {
146 BOOST_TEST_MESSAGE("Testing matrix exponential based on ode...");
147
148 // Reference results are taken from
149 // http://www.millersville.edu/~bikenaga/linear-algebra/matrix-exponential/matrix-exponential.html
150
151 Matrix m(3, 3);
152 m[0][0] = 5; m[0][1] =-6; m[0][2] =-6;
153 m[1][0] =-1; m[1][1] = 4; m[1][2] = 2;
154 m[2][0] = 3; m[2][1] =-6; m[2][2] =-4;
155
156 const Real tol = 1e-12;
157
158 for (Real t=0.01; t < 11; t+=t) {
159 const Matrix calculated = Expm(M: m, t, tol);
160
161 Matrix expected(3, 3);
162 expected[0][0] = -3*std::exp(x: t)+4*std::exp(x: 2*t);
163 expected[0][1] = 6*std::exp(x: t)-6*std::exp(x: 2*t);
164 expected[0][2] = 6*std::exp(x: t)-6*std::exp(x: 2*t);
165 expected[1][0] = std::exp(x: t)- std::exp(x: 2*t);
166 expected[1][1] = -2*std::exp(x: t)+3*std::exp(x: 2*t);
167 expected[1][2] = -2*std::exp(x: t)+2*std::exp(x: 2*t);
168 expected[2][0] = -3*std::exp(x: t)+3*std::exp(x: 2*t);
169 expected[2][1] = 6*std::exp(x: t)-6*std::exp(x: 2*t);
170 expected[2][2] = 6*std::exp(x: t)-5*std::exp(x: 2*t);
171
172 Matrix diff = calculated - expected;
173 Real relDiffNorm = frobenuiusNorm(m: diff)/frobenuiusNorm(m: expected);
174
175 if ( std::fabs(x: relDiffNorm) > 100*tol) {
176 BOOST_FAIL("Failed to reproduce expected matrix exponential."
177 << "\n rel. difference norm: " << relDiffNorm
178 << "\n tolerance : " << 100*tol);
179 }
180
181 const Matrix negativeTime = Expm(M: (-1)*m, t: -t, tol);
182 diff = negativeTime - expected;
183 relDiffNorm = frobenuiusNorm(m: diff)/frobenuiusNorm(m: expected);
184
185 if ( std::fabs(x: relDiffNorm) > 100*tol) {
186 BOOST_FAIL("Failed to reproduce expected matrix exponential."
187 << "\n rel. difference norm: " << relDiffNorm
188 << "\n tolerance : " << 100*tol);
189 }
190
191 }
192}
193
194void OdeTest::testMatrixExponentialOfZero() {
195 BOOST_TEST_MESSAGE("Testing matrix exponential of a zero matrix "
196 "based on ode...");
197
198 Matrix m(3, 3, 0.0);
199
200 constexpr double tol = 100*QL_EPSILON;
201 constexpr double t=1.0;
202 const Matrix calculated = Expm(M: m, t);
203
204 for (Size i=0; i < calculated.rows(); ++i) {
205 for (Size j=0; j < calculated.columns(); ++j) {
206 const Real kroneckerDelta = (i==j)? 1.0 : 0.0;
207 if (std::fabs(x: calculated[i][j] -kroneckerDelta) > tol) {
208 BOOST_FAIL("Failed to reproduce expected matrix exponential."
209 << "\n tolerance : " << tol);
210 }
211 }
212 }
213}
214
215test_suite* OdeTest::suite() {
216 auto* suite = BOOST_TEST_SUITE("ode tests");
217 suite->add(QUANTLIB_TEST_CASE(&OdeTest::testAdaptiveRungeKutta));
218 suite->add(QUANTLIB_TEST_CASE(&OdeTest::testMatrixExponential));
219 suite->add(QUANTLIB_TEST_CASE(&OdeTest::testMatrixExponentialOfZero));
220 return suite;
221}
222

source code of quantlib/test-suite/ode.cpp