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)); }
3.6.1 Finite Difference Schemes
We have written some simple code to implement the explicit Euler and Richardson extrapolation methods. We present it mainly for pedagogical reasons. First, we define global arrays that hold the discrete values of the three solutions (two Euler methods on two meshes and the second-order extrapolated scheme):
// Global arrays to hold (t, y) data for meshes of size dt and dt/2 // t and y arrays on mesh dt std::vector<double> tValues; std::vector<double> yValues; // t and y arrays on mesh dt/2 std::vector<double> t2Values; std::vector<double> y2Values; // 2nd order extrapolated values on mesh dt std::vector<double> yRichValues;
The code for the schemes is:
void CurveEuler() { // db/dt = F(b), b(0) given // b(n+1) = b(n) + dt * F(b(n)) tValues = CreateMesh(L, U, NT); yValues = std::vector<double>(NT+1); yValues[0] = A; auto riccatiEquation = CreateRiccati(); double yOld; double dt = (U - L) / static_cast<double>(NT); for (std::size_t j = 1; j < tValues.size(); ++j) { yOld = yValues[j-1]; yValues[j] = yOld + dt*(*riccatiEquation)(tValues[j-1], yOld); } } void CurveRichardson() { // O(dt) -> O(dt^2) // Exercise: remove duplicate code double dt = (U - L)/ double(NT); tValues = CreateMesh(L, U, NT); yValues = std::vector<double>(NT+1); yValues[0] = A; auto riccatiEquation = CreateRiccati(); double yOld; for (std::size_t j = 1; j < tValues.size(); ++j) { yOld = yValues[j-1]; yValues[j] = yOld + dt*(*riccatiEquation)(tValues[j - 1], yOld); } // dt/2 dt *= 0.5; t2Values = CreateMesh(L, U, 2*NT+1); y2Values = std::vector<double>(2*NT+1); y2Values[0] = A; for (std::size_t j = 1; j < t2Values.size(); ++j) { yOld = y2Values[j-1]; y2Values[j] = yOld + dt*(*riccatiEquation)(t2Values[j - 1], yOld); } yRichValues = std::vector<double>(NT+1); yRichValues[0] = A; for (std::size_t j = 1; j < yRichValues.size(); ++j) { yRichValues[j] = 2.0* y2Values[2*j] - yValues[j]; } }
Some test code is:
int main() { std::cout ≪ "How many steps (need many e.g. 300000?): "; std::cin ≫ NT; CurveEuler(); CurveRichardson(); std::cout ≪ "Value at T = " ≪ U ≪ " is " ≪ std::setprecision(7) ≪ yValues[yValues.size() - 1] ≪ ", press the ANY key " ≪ std::endl; std::cout ≪ "Richardson at T = " ≪ U ≪ " is " ≪RichValues[yRichValues.size() - 1] ≪ ", press the ANY key "; // Beautiful Excel output ExcelDriver xl; xl.MakeVisible(true); xl.CreateChart(tValues, yValues, std::string("Euler for Riccati")); xl.CreateChart (tValues, yRichValues, std::string("Richardson for Riccati")); return 0; }
What is the accuracy of these schemes, and how much computer effort (size of NT and the number of function calls) is needed in order to achieve a given accuracy? To this end, we take two cases whose exact solutions have been computed using Boost C++ odeint:
Case I: , .
Case II: , .
In general, the library needs approximately 100 function evaluations to reach this level of accuracy. The results for explicit Euler and the extrapolated scheme for cases I and II for NT = 300, 500, 1000 and 5000 are:
Case I: (9.2746e-6, 1.113e-5), (1.0015e-5, 1.118e-5), (1.059e-5, 1.200e-5), (1.083e-5, 1.1206e-5).
Case II: (0.0554, 0.0559), (0.0556, 0.0559), (0.0558, 0.0559), (0.05589, 0.0559).
For a more complete discussion of Boost odeint, see Duffy (2018).
We complete our discussion of finite difference methods for (3.10). The full problem is (3.11); (3.12) is a coupled system of equations, and it can be solved in different ways. The most effective approach at this stage is to use a production third-party ODE solver in C++.
3.7 MATRIX DIFFERENTIAL EQUATIONS
We have seen what scalar and vector ODEs are. Now we consider matrix ODEs.
The Boost odeint library can be used to solve matrix ODEs, for example:
where A, B and Y are
It can be shown that the solution of system (3.25) can be written as (see Brauer and Nohel (1969)):
(3.26)
A special case is when:
(3.27)
(3.28)
The conclusion is that we can now compute the exponential of a matrix by solving a matrix ODE. We now discuss this topic using C++. To this end, we examine the system:
(3.29)
where C is a given matrix.
We model this system by the following C++ class:
namespace ublas = boost::numeric::ublas; using value_type = double; using state_type = boost::numeric::ublas::matrix<value_type>; class MatrixOde { private: // dB/dt = A*B, B(0) = C; ublas::matrix<value_type> A_; ublas::matrix<value_type>
