funm

一般的な行列関数の計算

構文

F = funm(A,fun) F = funm(A,fun,options) F = funm(A,fun,options,p1,p2,...) [F,exitflag] = funm(...) [F,exitflag,output] = funm(...)

説明

F = funm(A,fun) は、ユーザー定義の関数 fun を正方行列の引数 A において評価します。F = fun(x,k) はベクトル x と整数 k を受け取り、x と同じサイズのベクトル f を返さなければなりません。ここで、f(i) は、x(i) で評価された関数 fun の k 番目の導関数です。fun で表される関数は、収束半径が無限大のテイラー級数をもつ必要があります。ただし、fun = @log の場合を除きます。これは、特殊なケースとして扱われます。

関数 funm を使用して、次の表に一覧表示された特殊な関数を行列 A で評価することもできます。

関数	関数を行列 A で評価するための構文
`exp`	`funm(A, @exp)`
`log`	`funm(A, @log)`
`sin`	`funm(A, @sin)`
`cos`	`funm(A, @cos)`
`sinh`	`funm(A, @sinh)`
`cosh`	`funm(A, @cosh)`

行列平方根の場合には、代わりに sqrtm(A) を使用します。行列指数の場合、expm(A) または funm(A, @exp) のいずれが正確かは、行列 A に依存します。

fun で表される関数は、収束半径が無限大のテイラー級数をもつ必要があります。例外は @log で、特殊なケースとして扱われます。必要な場合は、関数 fun に追加パラメーターを提供する方法が関数のパラメーター化に記載されています。

F = funm(A,fun,options) は、アルゴリズムのパラメーターを options 構造体の値に設定します。

次の表に、options のフィールドの一覧を示します。

フィールド	説明	値
`options.Display`	表示のレベル	`'off'` (既定の設定)、`'on'`、`'verbose'`
`options.TolBlk`	Schur 型のブロックに対する許容誤差	正のスカラー。既定値は `0.1`です。
`options.TolTay`	対角ブロックのテイラー級数を評価するための終了許容誤差	正のスカラー。既定値は `eps` です。
`options.MaxTerms`	テイラー級数の最大項数	正の整数。既定値は `250`です。
`options.MaxSqrt`	アルゴリズムの計算において、逆スケーリングと二乗法で計算された平方根の最大数	正の整数。既定値は `100`です。
`options.Ord`	Schur 型 `T` の順序を指定します。	長さが `length(A)` のベクトル。`options.Ord(i)` は `T(i,i)` が配置されるブロックのインデックスです。既定値は `[]` です。

F = funm(A,fun,options,p1,p2,...) は、余分の入力 p1,p2,... を関数に渡します。

[F,exitflag] = funm(...) は、関数 funm の終了条件を記述するスカラー exitflag を返します。引数 exitflag は次の値を取ります。

0 — アルゴリズムは成功しました。
1 — 1 つ以上のテイラー級数が収束しなかった、または、対数の場合は、必要な行列の平方根が多すぎたことを示します。ただし、その場合にも F の値は正しく計算されている可能性があります。

[F,exitflag,output] = funm(...) は、以下のフィールドをもつ構造体 output を出力します。

フィールド	説明
`output.terms`	`output.terms(i)` が i 番目のブロックの評価に使用する際のテイラー級数の項数となるベクトル、または対数の場合には、2 より大きい次元の行列の平方根の数
`output.ind`	並べ替えられた Schur 係数 `T` の `(i,j)` ブロックが `T(output.ind{i}, output.ind{j})` となる cell 配列
`output.ord`	Schur 型が関数 `ordschur` に渡されるときの順序
`output.T`	並べ替えた Schur 型

Schur が対角行列の場合、output = struct('terms',ones(n,1),'ind',{1:n}) となります。

例

例 1

次のコマンドは、3 行 3 列の魔方陣行列の正弦値を計算します。

F=funm(magic(3), @sin)

F =

   -0.3850    1.0191    0.0162
    0.6179    0.2168   -0.1844
    0.4173   -0.5856    0.8185

例 2

ステートメント

S = funm(X,@sin);
C = funm(X,@cos);

は、丸め誤差の範囲内で以下の場合と同じ結果を出力します。

E = expm(i*X);
C = real(E);
S = imag(E);

いずれの場合も結果は S*S+C*C = I を満たします。ここで、I = eye(size(X)) です。

例 3

A における関数 exp(x) + cos(x) を 1 回の関数 funm の呼び出しで計算する方法は以下のとおりです。

F = funm(A,@fun_expcos)

ここで、fun_expcos は次の関数です。

function f = fun_expcos(x, k)
% Return kth derivative of exp + cos at X.
        g = mod(ceil(k/2),2);
        if mod(k,2)
           f = exp(x) + sin(x)*(-1)^g;
        else
           f = exp(x) + cos(x)*(-1)^g;
        end

アルゴリズム

関数 funm のアルゴリズムは[1]に記述されています。

参照

[1] Davies, P. I. and N. J. Higham, “A Schur-Parlett algorithm for computing matrix functions,” SIAM J. Matrix Anal. Appl., Vol. 25, Number 2, pp. 464-485, 2003.

[2] Golub, G. H. and C. F. Van Loan, Matrix Computation, Third Edition, Johns Hopkins University Press, 1996, p. 384.

[3] Moler, C. B. and C. F. Van Loan, “Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later” SIAM Review 20, Vol. 45, Number 1, pp. 1-47, 2003.