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

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:

(3.25)

where A, B and Y are matrices.

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> C_; public: MatrixOde(const ublas::matrix<value_type>& A, const ublas::matrix<value_type>& IC) : A_(A), C_(IC) {} void operator()(const state_type &x , state_type &dxdt, double t ) const { for( std::size_t i=0 ; i < x.size1();++i ) { for( std::size_t j=0 ; j < x.size2(); ++j ) { dxdt(i, j) = 0.0; for (std::size_t k = 0; k < x.size2(); ++k) { dxdt(i, j) += A_(i,k)*x(k,j); } } } } };

There are many dubious ways to compute the exponential of a matrix, see Moler and Van Loan (2003), one of which involves the application of ODE solvers. Other methods include:

 S1: Series methods (for example, truncating the infinite Taylor series representation for the exponential).

 S2: Padé rational approximant. This entails approximating the exponential by a special kind of rational function.

 S3: Polynomial methods using the Cayley–Hamilton method.

 S4: Inverse Laplace transform.

 S5: Matrix decomposition methods.

 S6: Splitting methods.