1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2012 Peter Caspers
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 adaptiverungekutta.hpp
21 \brief Runge-Kutta ODE integration
22
23 Runge Kutta method with adaptive stepsize as described in
24 Numerical Recipes in C, Chapter 16.2
25*/
26
27#ifndef quantlib_adaptive_runge_kutta_hpp
28#define quantlib_adaptive_runge_kutta_hpp
29
30#include <ql/types.hpp>
31#include <ql/errors.hpp>
32#include <ql/functional.hpp>
33#include <vector>
34#include <cmath>
35#include <complex>
36
37namespace QuantLib {
38
39 template <class T = Real>
40 class AdaptiveRungeKutta {
41 public:
42 typedef ext::function<std::vector<T>(const Real, const std::vector<T>&)> OdeFct;
43 typedef ext::function<T(const Real, const T)> OdeFct1d;
44
45 /*! The class is constructed with the following inputs:
46 - eps prescribed error for the solution
47 - h1 start step size
48 - hmin smallest step size allowed
49 */
50
51 AdaptiveRungeKutta(const Real eps = 1.0e-6, const Real h1 = 1.0e-4, const Real hmin = 0.0)
52 : eps_(eps), h1_(h1), hmin_(hmin), b31(3.0 / 40.0), b32(9.0 / 40.0), b51(-11.0 / 54.0),
53 b53(-70.0 / 27.0), b54(35.0 / 27.0), b61(1631.0 / 55296.0), b62(175.0 / 512.0),
54 b63(575.0 / 13824.0), b64(44275.0 / 110592.0), b65(253.0 / 4096.0), c1(37.0 / 378.0),
55 c3(250.0 / 621.0), c4(125.0 / 594.0), c6(512.0 / 1771.0), dc1(c1 - 2825.0 / 27648.0),
56 dc3(c3 - 18575.0 / 48384.0), dc4(c4 - 13525.0 / 55296.0), dc5(-277.0 / 14336.0),
57 dc6(c6 - 0.25) {}
58
59 /*! Integrate the ode from \f$ x1 \f$ to \f$ x2 \f$ with
60 initial value condition \f$ f(x1)=y1 \f$.
61
62 The ode is given by a function \f$ F: R \times K^n
63 \rightarrow K^n \f$ as \f$ f'(x) = F(x,f(x)) \f$, $K=R,
64 C$ */
65 std::vector<T> operator()(const OdeFct& ode, const std::vector<T>& y1, Real x1, Real x2);
66 T operator()(const OdeFct1d& ode, T y1, Real x1, Real x2);
67
68 private:
69 void rkqs(std::vector<T>& y,
70 const std::vector<T>& dydx,
71 Real& x,
72 Real htry,
73 Real eps,
74 const std::vector<Real>& yScale,
75 Real& hdid,
76 Real& hnext,
77 const OdeFct& derivs);
78 void rkck(const std::vector<T>& y,
79 const std::vector<T>& dydx,
80 Real x,
81 Real h,
82 std::vector<T>& yout,
83 std::vector<T>& yerr,
84 const OdeFct& derivs);
85
86 const std::vector<T> yStart_;
87 const Real eps_, h1_, hmin_;
88 const Real a2 = 0.2, a3 = 0.3, a4 = 0.6, a5 = 1.0, a6 = 0.875, b21 = 0.2, b31, b32,
89 b41 = 0.3, b42 = -0.9, b43 = 1.2, b51, b52 = 2.5, b53, b54, b61, b62, b63, b64,
90 b65, c1, c3, c4, c6, dc1, dc3, dc4, dc5, dc6;
91 const double ADAPTIVERK_MAXSTP = 10000, ADAPTIVERK_TINY = 1.0E-30, ADAPTIVERK_SAFETY = 0.9,
92 ADAPTIVERK_PGROW = -0.2, ADAPTIVERK_PSHRINK = -0.25,
93 ADAPTIVERK_ERRCON = 1.89E-4;
94 };
95
96
97
98 template<class T>
99 std::vector<T> AdaptiveRungeKutta<T>::operator()(const OdeFct& ode,
100 const std::vector<T>& y1,
101 const Real x1,
102 const Real x2) {
103 Size n = y1.size();
104 std::vector<T> y(y1);
105 std::vector<Real> yScale(n);
106 Real x = x1;
107 Real h = h1_* (x1<=x2 ? 1 : -1);
108 Real hnext,hdid;
109
110 for (Size nstp=1; nstp<=ADAPTIVERK_MAXSTP; nstp++) {
111 std::vector<T> dydx=ode(x,y);
112 for (Size i=0;i<n;i++)
113 yScale[i] = std::abs(y[i])+std::abs(dydx[i]*h)+ADAPTIVERK_TINY;
114 if ((x+h-x2)*(x+h-x1) > 0.0)
115 h=x2-x;
116 rkqs(y,dydx,x,htry: h,eps: eps_,yScale,hdid,hnext,derivs: ode);
117
118 if ((x-x2)*(x2-x1) >= 0.0)
119 return y;
120
121 if (std::fabs(x: hnext) <= hmin_)
122 QL_FAIL("Step size (" << hnext << ") too small ("
123 << hmin_ << " min) in AdaptiveRungeKutta");
124 h=hnext;
125 }
126 QL_FAIL("Too many steps (" << ADAPTIVERK_MAXSTP
127 << ") in AdaptiveRungeKutta");
128 }
129
130 namespace detail {
131
132 template <class T>
133 struct OdeFctWrapper {
134 typedef typename AdaptiveRungeKutta<T>::OdeFct1d OdeFct1d;
135 explicit OdeFctWrapper(const OdeFct1d& ode1d)
136 : ode1d_(ode1d) {}
137 std::vector<T> operator()(const Real x, const std::vector<T>& y) {
138 std::vector<T> res(1,ode1d_(x,y[0]));
139 return res;
140 }
141 const OdeFct1d& ode1d_;
142 };
143
144 }
145
146 template<class T>
147 T AdaptiveRungeKutta<T>::operator()(const OdeFct1d& ode,
148 const T y1,
149 const Real x1,
150 const Real x2) {
151 return operator()(detail::OdeFctWrapper<T>(ode),
152 std::vector<T>(1,y1),x1,x2)[0];
153 }
154
155 template<class T>
156 void AdaptiveRungeKutta<T>::rkqs(std::vector<T>& y,
157 const std::vector<T>& dydx,
158 Real& x,
159 const Real htry,
160 const Real eps,
161 const std::vector<Real>& yScale,
162 Real& hdid,
163 Real& hnext,
164 const OdeFct& derivs) {
165 Size n=y.size();
166 Real errmax,xnew;
167 std::vector<T> yerr(n),ytemp(n);
168
169 Real h=htry;
170
171 for(;;) {
172 rkck(y,dydx,x,h,yout&: ytemp,yerr,derivs);
173 errmax=0.0;
174 for (Size i=0;i<n;i++)
175 errmax=std::max(errmax,std::abs(yerr[i]/yScale[i]));
176 errmax/=eps;
177 if (errmax>1.0) {
178 Real htemp1 = ADAPTIVERK_SAFETY*h*std::pow(x: errmax,y: ADAPTIVERK_PSHRINK);
179 Real htemp2 = h / 10;
180 // These would be std::min and std::max, of course,
181 // but VC++14 had problems inlining them and caused
182 // the wrong results to be calculated. The problem
183 // seems to be fixed in update 3, but let's keep this
184 // implementation for compatibility.
185 Real max_positive = htemp1 > htemp2 ? htemp1 : htemp2;
186 Real max_negative = htemp1 < htemp2 ? htemp1 : htemp2;
187 h = ((h >= 0.0) ? max_positive : max_negative);
188 xnew=x+h;
189 if (xnew==x)
190 QL_FAIL("Stepsize underflow (" << h << " at x = " << x
191 << ") in AdaptiveRungeKutta::rkqs");
192 continue;
193 } else {
194 if (errmax>ADAPTIVERK_ERRCON)
195 hnext=ADAPTIVERK_SAFETY*h*std::pow(x: errmax,y: ADAPTIVERK_PGROW);
196 else
197 hnext=5.0*h;
198 x+=(hdid=h);
199 for (Size i=0;i<n;i++)
200 y[i]=ytemp[i];
201 break;
202 }
203 }
204 }
205
206 template <class T>
207 void AdaptiveRungeKutta<T>::rkck(const std::vector<T>& y,
208 const std::vector<T>& dydx,
209 Real x,
210 const Real h,
211 std::vector<T>& yout,
212 std::vector<T> &yerr,
213 const OdeFct& derivs) {
214
215 Size n=y.size();
216 std::vector<T> ak2(n),ak3(n),ak4(n),ak5(n),ak6(n),ytemp(n);
217
218 // first step
219 for (Size i=0;i<n;i++)
220 ytemp[i]=y[i]+b21*h*dydx[i];
221
222 // second step
223 ak2=derivs(x+a2*h,ytemp);
224 for (Size i=0;i<n;i++)
225 ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
226
227 // third step
228 ak3=derivs(x+a3*h,ytemp);
229 for (Size i=0;i<n;i++)
230 ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
231
232 // fourth step
233 ak4=derivs(x+a4*h,ytemp);
234 for (Size i=0;i<n;i++)
235 ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
236
237 // fifth step
238 ak5=derivs(x+a5*h,ytemp);
239 for (Size i=0;i<n;i++)
240 ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
241
242 // sixth step
243 ak6=derivs(x+a6*h,ytemp);
244 for (Size i=0;i<n;i++) {
245 yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
246 yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
247 }
248 }
249
250}
251
252#endif
253

source code of quantlib/ql/math/ode/adaptiverungekutta.hpp