Читать книгу Numerical Methods in Computational Finance - Daniel J. Duffy - Страница 62

We return to the non-linear IVP (3.10) that we solve numerically using the explicit Euler method, Richardson extrapolation, and several robust ODE solvers from the Boost C++ library odeint. We expect the explicit Euler scheme to perform badly because it is explicit and Equation (3.10) is non-linear. We describe the simple software design for the Euler and extrapolation schemes (the main objective is to show how to map numerical algorithms to code). First, we create a class that models the Riccati Equation (3.10):

using value_type = double; using FunctionType = std::function<value_type (value_type)>; using FunctionType2 = std::function<value_type (value_type x, value_type y)>; class GeneralisedRiccati { private: FunctionType p, q, r; FunctionType2 n; public: GeneralisedRiccati(const FunctionType& P, const FunctionType& Q, const FunctionType& R, const FunctionType2& N) : p(P), q(Q), r(R), n(N) { } double operator() (double x, double y) { // Function object (functor) return p(x) + q(x)*y + r(x)*y*y + n(x, y); } };

We can instantiate this class by calling its constructor, but we prefer to use a factory method (this is a design pattern) for added flexibility. Notice that we also use C++ smart pointers in order to avoid possible memory issues):

std::shared_ptr<GeneralisedRiccati> CreateRiccati() { // Factory method, specific case const double P = 0.0; const double Q = -10.8; const double R = 0.5 * 2.44 * 2.44; const double eta = 0.96; const double lambda = 0.13;double A = 0.42; // Generalised Riccati auto p = [=](double x) {return P;};auto q = [=](double x) {return Q;}; auto r = [=](double x) {return R;}; auto n = [=](double x, double y) { return (lambda*y) / (eta - y); }; return std::shared_ptr<GeneralisedRiccati> (new GeneralisedRiccati (p, q, r, n)); }