27#ifndef quantlib_adaptive_runge_kutta_hpp
28#define quantlib_adaptive_runge_kutta_hpp
39 template <
class T = Real>
42 typedef ext::function<std::vector<T>(
const Real,
const std::vector<T>&)>
OdeFct;
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),
69 void rkqs(std::vector<T>&
y,
70 const std::vector<T>& dydx,
74 const std::vector<Real>& yScale,
78 void rkck(
const std::vector<T>&
y,
79 const std::vector<T>& dydx,
89 b41 = 0.3,
b42 = -0.9,
b43 = 1.2,
b51,
b52 = 2.5,
b53,
b54,
b61,
b62,
b63,
b64,
100 const std::vector<T>& y1,
104 std::vector<T>
y(y1);
105 std::vector<Real> yScale(
n);
107 Real h = h1_* (x1<=x2 ? 1 : -1);
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)
116 rkqs(
y,dydx,x,h,eps_,yScale,hdid,hnext,ode);
118 if ((x-x2)*(x2-x1) >= 0.0)
121 if (std::fabs(hnext) <= hmin_)
122 QL_FAIL(
"Step size (" << hnext <<
") too small ("
123 << hmin_ <<
" min) in AdaptiveRungeKutta");
126 QL_FAIL(
"Too many steps (" << ADAPTIVERK_MAXSTP
127 <<
") in AdaptiveRungeKutta");
138 std::vector<T> res(1,
ode1d_(x,
y[0]));
152 std::vector<T>(1,y1),x1,x2)[0];
157 const std::vector<T>& dydx,
161 const std::vector<Real>& yScale,
167 std::vector<T> yerr(
n),ytemp(
n);
172 rkck(
y,dydx,x,h,ytemp,yerr,derivs);
174 for (
Size i=0;i<
n;i++)
175 errmax=std::max(errmax,std::abs(yerr[i]/yScale[i]));
178 Real htemp1 = ADAPTIVERK_SAFETY*h*std::pow(errmax,ADAPTIVERK_PSHRINK);
179 Real htemp2 = h / 10;
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);
190 QL_FAIL(
"Stepsize underflow (" << h <<
" at x = " << x
191 <<
") in AdaptiveRungeKutta::rkqs");
194 if (errmax>ADAPTIVERK_ERRCON)
195 hnext=ADAPTIVERK_SAFETY*h*std::pow(errmax,ADAPTIVERK_PGROW);
199 for (
Size i=0;i<
n;i++)
208 const std::vector<T>& dydx,
211 std::vector<T>& yout,
212 std::vector<T> &yerr,
216 std::vector<T> ak2(
n),ak3(
n),ak4(
n),ak5(
n),ak6(
n),ytemp(
n);
219 for (
Size i=0;i<
n;i++)
220 ytemp[i]=
y[i]+b21*h*dydx[i];
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]);
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]);
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]);
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]);
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]);
void rkqs(std::vector< T > &y, const std::vector< T > &dydx, Real &x, Real htry, Real eps, const std::vector< Real > &yScale, Real &hdid, Real &hnext, const OdeFct &derivs)
const double ADAPTIVERK_PSHRINK
const double ADAPTIVERK_PGROW
const double ADAPTIVERK_SAFETY
ext::function< std::vector< T >(const Real, const std::vector< T > &)> OdeFct
const double ADAPTIVERK_ERRCON
const double ADAPTIVERK_MAXSTP
void rkck(const std::vector< T > &y, const std::vector< T > &dydx, Real x, Real h, std::vector< T > &yout, std::vector< T > &yerr, const OdeFct &derivs)
const std::vector< T > yStart_
const double ADAPTIVERK_TINY
ext::function< T(const Real, const T)> OdeFct1d
std::vector< T > operator()(const OdeFct &ode, const std::vector< T > &y1, Real x1, Real x2)
AdaptiveRungeKutta(const Real eps=1.0e-6, const Real h1=1.0e-4, const Real hmin=0.0)
Classes and functions for error handling.
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Maps function, bind and cref to either the boost or std implementation.
std::size_t Size
size of a container
AdaptiveRungeKutta< T >::OdeFct1d OdeFct1d
OdeFctWrapper(const OdeFct1d &ode1d)
std::vector< T > operator()(const Real x, const std::vector< T > &y)