Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

シンボリック行列演算

この例では、Symbolic Math Toolbox™ を使って簡単な行列計算を実行する方法を示します。

一般的なテスト行列である 5 行 5 列のヒルベルト行列を生成します。

H = sym(hilb(5)) 
H = 

(1121314151213141516131415161714151617181516171819)

この行列式は非常に小さいです。

d = det(H) 
d = 

1266716800000

逆行列の要素は整数です。

X = inv(H) 
X = 

(25-3001050-1400630-3004800-1890026880-126001050-1890079380-11760056700-140026880-117600179200-88200630-1260056700-8820044100)

この逆行列が正しいことを確認します。

I = X*H
I = 

(1000001000001000001000001)

特性多項式を求めます。

syms x; p = charpoly(H,x) 
p = 

x5-563x4315+735781x32116800-852401x2222264000+61501x53343360000-1266716800000

この特性多項式の因数分解を試みます。

factor(p) 
ans = 

(1266716800000266716800000x5-476703360000x4+92708406000x3-1022881200x2+307505x-1)

この結果は、この特性多項式を有理数上で因数分解できないことを示しています。

固有値を 50 桁の精度で近似計算します。

digits(50) 
e = eig(vpa(H)) 
e = 

(1.56705069109823079553301100552072463394931525223340.208534218611013335905002510068820055038582022603430.011407491623419806559451458866589345042348430526640.000305898040151191726879497840692722825656145149092470.0000032879287721718629571150047605447313997367890230746)

自由変数 t を含む一般化ヒルベルト行列を作成します。

t = sym('t'); 
[I,J] = meshgrid(1:5); 
H = 1./(I+J-t)
H = 

(-1t-2-1t-3σ5σ3σ1-1t-3σ5σ3σ1σ2σ5σ3σ1σ2σ4σ3σ1σ2σ4-1t-9σ1σ2σ4-1t-9-1t-10)where  σ1=-1t-6  σ2=-1t-7  σ3=-1t-5  σ4=-1t-8  σ5=-1t-4

t=1 を代入すると、元のヒルベルト行列が得られます。

subs(H,t,1) 
ans = 

(1121314151213141516131415161714151617181516171819)

行列式の逆数は、t の多項式です。

d = 1/det(H) 
d = 

-t-2t-32t-43t-54t-65t-74t-83t-92t-1082944

d = expand(d)
d = 

-t2582944+25t2413824-5375t2341472+40825t226912-15940015t2182944+21896665t204608-240519875t192592+1268467075t18864-1588946776255t1782944+2885896606895t1613824-79493630114675t1541472+34372691161375t142304-8194259295156385t1382944+7707965729450845t1213824-55608098247105175t1120736+37909434298793825t103456-197019820623693025t95184+10640296363350955t896-38821472549340925t7144+12958201048605475t624-1748754621252377t52+1115685328012530t4-1078920141906600t3+742618453752000t2-323874210240000t+67212633600000

逆行列の要素も t の多項式です。

X = inv(H) 
X = 

(-t-2t-3t-4t-5t-6σ4576t-3t-4t-5t-6t-7t4-17t3+104t2-268t+240144-t-4t-5t-6t-7t-8t4-16t3+91t2-216t+18096t-5t-6t-7t-8t-9t4-15t3+80t2-180t+144144-t-6t-7t-8t-9t-10t4-14t3+71t2-154t+120576t-2t-3t-4t-5t-6σ3144-t-3t-4t-5t-6t-7t4-21t3+161t2-531t+63036t-4t-5t-6t-7t-8t4-20t3+145t2-450t+50424-t-5t-6t-7t-8t-9t4-19t3+131t2-389t+42036t-6t-7t-8t-9t-10σ4144-t-2t-3t-4t-5t-6σ296t-3t-4t-5t-6t-7t4-25t3+230t2-920t+134424-t-4t-5t-6t-7t-8t4-24t3+211t2-804t+112016t-5t-6t-7t-8t-9t4-23t3+194t2-712t+96024-t-6t-7t-8t-9t-10σ396t-2t-3t-4t-5t-6σ1144-t-3t-4t-5t-6t-7t4-29t3+311t2-1459t+252036t-4t-5t-6t-7t-8t4-28t3+289t2-1302t+216024-t-5t-6t-7t-8t-9t4-27t3+269t2-1173t+189036t-6t-7t-8t-9t-10σ2144-t-2t-3t-4t-5t-6t4-34t3+431t2-2414t+5040576t-3t-4t-5t-6t-7t4-33t3+404t2-2172t+4320144-t-4t-5t-6t-7t-8t4-32t3+379t2-1968t+378096t-5t-6t-7t-8t-9t4-31t3+356t2-1796t+3360144-t-6t-7t-8t-9t-10σ1576)where  σ1=t4-30t3+335t2-1650t+3024  σ2=t4-26t3+251t2-1066t+1680  σ3=t4-22t3+179t2-638t+840  σ4=t4-18t3+119t2-342t+360

t=1 を代入すると、ヒルベルト行列の逆行列が生成されます。

X = subs(X,t,'1') 
X = 

(25-3001050-1400630-3004800-1890026880-126001050-1890079380-11760056700-140026880-117600179200-88200630-1260056700-8820044100)

X = double(X) 
X = 5×5

          25        -300        1050       -1400         630
        -300        4800      -18900       26880      -12600
        1050      -18900       79380     -117600       56700
       -1400       26880     -117600      179200      -88200
         630      -12600       56700      -88200       44100

別の例を調べます。

A = sym(gallery(5)) 
A = 

(-911-2163-25270-69141-4211684-575575-11493451-138013891-38917782-23345933651024-10242048-614424572)

これは "べき零行列" です。これを 5 乗すると、ゼロ行列になります。

A^5 
ans = 

(0000000000000000000000000)

これはべき零行列であるため、その特性多項式は非常に単純です。

p = charpoly(A,'lambda') 
p = λ5

ここで、この固有値の計算は暗算ができるぐらい単純です。固有値は、方程式 lambda^5 = 0 の零点です。

シンボリック計算によって、正確な固有値を求めることができます。

lambda = eig(A) 
lambda = 

(00000)

数値計算は丸め誤差を伴い、lambda^5 = eps*norm(A) のような方程式の零点を求めます。このため、計算された固有値は概ね lambda = (eps*norm(A))^(1/5) となります。ここで、Symbolic Math Toolbox によって 16 桁の精度で、固有値の計算をします。すべての解が 0 になるべきですが、計算結果を確認すると、そうではないことが分かります。

digits(16) 
lambda = eig(vpa(A)) 
lambda = 

(0.00056174484863958470.0001737348850386136-0.0005342985684139864i0.0001737348850386136+0.0005342985684139864i-0.000454607309358406+0.0003303865815979566i-0.000454607309358406-0.0003303865815979566i)

これは、"欠陥" 行列でもあります。この行列は、対数行列ではありません。この行列のジョルダン正準形は非対角です。

J = jordan(A) 
J = 

(0100000100000100000100000)

行列の指数関数 expm(t*A) は通常、固有値 exp(lambda(i)*t) を含むスカラーの指数関数によって表されます。しかし、この行列の場合、expm(t*A) の要素はすべて t の多項式です。

t = sym('t'); 
E = simplify(expm(t*A)) 
E = 

(-2t33+11t22-9t+1t4t2-27t+333-t20t2-117t+1266t32t2-174t+1893-t85t2-464t+5042-t7t3-81t2+230t-14027t4-67t3+301t22-69t+1-t35t3-293t2+598t-2822t112t3-876t2+1799t-8422-t1785t3-14012t2+28776t-134728t142t3-1710t2+5151t-34506-t142t3-1426t2+3438t-17253355t43-3139t33+4585t22-1149t+1-t1136t3-9420t2+20625t-103533t18105t3-150646t2+329952t-16561212-t973t3-11675t2+35022t-233466t1946t3-19458t2+46695t-233466-t4865t3-42807t2+93390t-4669267784t43-64210t33+93391t22-23345t+1-t248115t3-2053748t2+4482036t-224076024-128tt3-12t2+36t-243256tt3-10t2+24t-123-128t5t3-44t2+96t-483512t4t3-33t2+72t-363-2720t4+67552t33-49144t2+24572t+1)

ところで、関数 "exp" は要素ごとの指数を計算します。

X = exp(t*A) 
X = 

(e-9te11te-21te63te-252te70te-69te141te-421te1684te-575te575te-1149te3451te-13801te3891te-3891te7782te-23345te93365te1024te-1024te2048te-6144te24572t)