EigenのMatrix function機能

はじめに

EigenにはMatrix functions moduleがあります。これを使うことで、行列のcosやexpを計算することができます。

例(コサイン)

そもそも行列のcosって何でしょう。「matrix function cos」で検索すると、こちらのサイトが見つかりました。
What is the cosine of a matrix? | John D. Cook
このページによると、テイラー展開を行うことで計算をするようです。EigenとMATLABで実装した例を以下に載せます。

C++

MatrixXd A(2, 2);
A << 1, 2,
	 3, 4;
std::cout << A.cos() << std::endl;
// 結果
//  0.855423 -0.110876
// -0.166315  0.689109

MATLAB

A = [1 2; 3 4]
funm(A, @cos)

I = eye(2);
I - A^2/factorial(2) ...
  + A^4/factorial(4) ...
  - A^6/factorial(6) ...
  + A^8/factorial(8) ...
  - A^10/factorial(10) ...
  + A^12/factorial(12) ...
  - A^14/factorial(14) ...

MATLABの出力

A =
     1     2
     3     4

ans =
    0.8554   -0.1109
   -0.1663    0.6891

ans =
    0.8504   -0.1182
   -0.1773    0.6731

例(exp)

もう一つの例として行列のexpを計算してみます。

これは現代制御の分野では遷移行列といわれ有名なようです。（MATLAB/Simulinkによる現代制御入門）。せっかくなのでこの本の例題を再現してみます。
現代制御では以下の式で線形システムを表現できます。

バネとダンパにつながった台車を動かすときの応答を、バネの力など考慮することで、

と具体化しています。入力u(t)=0、つまりゼロ入力応答は

より

と解が求まります。ここで遷移行列が現れました。時刻tから位置y(t)が算出できます。
MATLAB

A = [0 1 
    -10 -2];
C = [1 0];

x0 = [1
    0];
t = 1;
yt = C*expm(A*t)*x0

MATLABの実行結果

yt =
   -0.3469

C++

#include <iostream>
#include <unsupported/Eigen/MatrixFunctions>
using namespace Eigen;

std::complex<double> expfn(std::complex<double> x, int)
{
  return std::exp(x);
}

int main()
{
  const double t = 1.0;
  MatrixXd A(2, 2);
  MatrixXd c(1, 2);
  VectorXd x0(2);

  A << 0, 1, -10, -2;
  c << 1, 0;
  x0 << 1, 0;

  std::cout << c * (A*t).matrixFunction(expfn) * x0;
 // 結果：-0.346893
}

はじめ1mだけ台車をずらしておいて、手を離すとバネに押されて戻るイメージです。