ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

行列の指数

この例では、行列の指数を計算する 19 通りの方法のうち 3 つを挙げています。

行列の指数の計算の背景については、SIAM Review 45, 3-49, 2003 の "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later" を参照してください。

Pseudospectra Gateway も推奨します。Web サイトは、以下のとおりです。

http://web.comlab.ox.ac.uk/projects/pseudospectra/

行列 A で開始

A = [0 1 2; 0.5 0 1; 2 1 0]
Asave = A;
A =

         0    1.0000    2.0000
    0.5000         0    1.0000
    2.0000    1.0000         0

スケーリングと正方化

expmdemo1 は、Golub and Van Loan, Matrix Computations, 3rd edition のアルゴリズム 11.3.1 を実装したものです。

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
   c = c * (q-k+1) / (k*(2*q-k+1));
   X = A*X;
   cX = c*X;
   E = E + cX;
   if p
     D = D + cX;
   else
     D = D - cX;
   end
   p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
for k = 1:s
    E = E*E;
end

E1 = E
E1 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

テイラー級数による行列の指数

expmdemo2 は、行列の指数に対する古典的な定義を使用します。実用的な数値的方法として、norm(A) が大きすぎる場合、これは遅い上に不正確です。

A = Asave;

% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
   E = E + F;
   F = A*F/k;
   k = k+1;
end

E2 = E
E2 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

固有値と固有ベクトルによる行列の指数

expmdemo3 は、行列が固有ベクトルの完全なセットをもつものと仮定します。実際の数値的方法として、精度は固有ベクトル行列の条件で決まります。

A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E
E3 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

結果の比較

この行列に対して、同じように実行されます。

E = expm(Asave);
err1 = E - E1
err2 = E - E2
err3 = E - E3
err1 =

   1.0e-14 *

    0.3553    0.1776    0.0888
    0.0888    0.1332   -0.0444
         0         0   -0.2665


err2 =

   1.0e-14 *

         0         0   -0.1776
   -0.0444         0   -0.0888
    0.1776         0    0.0888


err3 =

   1.0e-13 *

   -0.0799   -0.0444   -0.0622
   -0.0622   -0.0488   -0.0933
   -0.0711   -0.0533   -0.1066

テイラー級数の失敗

これらが 0 になる前にテイラー級数の項が非常に大きくなる行列があります。したがって、expmdemo2 は失敗します。

A = [-147 72; -192 93];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)
E1 =

   -0.0996    0.0747
   -0.1991    0.1494


E2 =

   1.0e+06 *

   -1.1985   -0.5908
   -2.7438   -2.0442


E3 =

   -0.0996    0.0747
   -0.1991    0.1494

固有値とベクトルの失敗

固有ベクトルの完全なセットをもたない行列があります。したがって、expmdemo3 は失敗します。

A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)
E1 =

    0.3679    0.3679
         0    0.3679


E2 =

    0.3679    0.3679
         0    0.3679


E3 =

    0.3679         0
         0    0.3679

この情報は役に立ちましたか?