Main Content

このページは前リリースの情報です。該当の英語のページはこのリリースで削除されています。

シンボリック行列演算

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

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

H = sym(hilb(5)) 
H = 

(1121314151213141516131415161714151617181516171819)[sym(1), sym(1/2), sym(1/3), sym(1/4), sym(1/5); sym(1/2), sym(1/3), sym(1/4), sym(1/5), sym(1/6); sym(1/3), sym(1/4), sym(1/5), sym(1/6), sym(1/7); sym(1/4), sym(1/5), sym(1/6), sym(1/7), sym(1/8); sym(1/5), sym(1/6), sym(1/7), sym(1/8), sym(1/9)]

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

d = det(H) 
d = 

1266716800000sym('1/266716800000')

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

X = inv(H) 
X = 

(25-3001050-1400630-3004800-1890026880-126001050-1890079380-11760056700-140026880-117600179200-88200630-1260056700-8820044100)[sym(25), -sym(300), sym(1050), -sym(1400), sym(630); -sym(300), sym(4800), -sym(18900), sym(26880), -sym(12600); sym(1050), -sym(18900), sym(79380), -sym(117600), sym(56700); -sym(1400), sym(26880), -sym(117600), sym(179200), -sym(88200); sym(630), -sym(12600), sym(56700), -sym(88200), sym(44100)]

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

I = X*H
I = 

(1000001000001000001000001)[sym(1), sym(0), sym(0), sym(0), sym(0); sym(0), sym(1), sym(0), sym(0), sym(0); sym(0), sym(0), sym(1), sym(0), sym(0); sym(0), sym(0), sym(0), sym(1), sym(0); sym(0), sym(0), sym(0), sym(0), sym(1)]

特性多項式を求めます。

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

x5-563x4315+735781x32116800-852401x2222264000+61501x53343360000-1266716800000x^5 - (563*x^4)/315 + (735781*x^3)/2116800 - (852401*x^2)/222264000 + (61501*x)/sym('53343360000') - sym('1/266716800000')

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

factor(p) 
ans = 

(1266716800000266716800000x5-476703360000x4+92708406000x3-1022881200x2+307505x-1)[sym('1/266716800000'), sym('266716800000')*x^5 - sym('476703360000')*x^4 + sym('92708406000')*x^3 - 1022881200*x^2 + 307505*x - 1]

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

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

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

(1.56705069109823079553301100552072463394931525223340.208534218611013335905002510068820055038582022603430.011407491623419806559451458866589345042348430526640.000305898040151191726879497840692722825656145149092470.0000032879287721718629571150047605447313997367890230746)[vpa('1.5670506910982307955330110055207246339493152522334'); vpa('0.20853421861101333590500251006882005503858202260343'); vpa('0.01140749162341980655945145886658934504234843052664'); vpa('0.00030589804015119172687949784069272282565614514909247'); vpa('0.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[-1/(t - 2), -1/(t - 3), -1/(t - 4), -1/(t - 5), -1/(t - 6); -1/(t - 3), -1/(t - 4), -1/(t - 5), -1/(t - 6), -1/(t - 7); -1/(t - 4), -1/(t - 5), -1/(t - 6), -1/(t - 7), -1/(t - 8); -1/(t - 5), -1/(t - 6), -1/(t - 7), -1/(t - 8), -1/(t - 9); -1/(t - 6), -1/(t - 7), -1/(t - 8), -1/(t - 9), -1/(t - 10)]

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

subs(H,t,1) 
ans = 

(1121314151213141516131415161714151617181516171819)[sym(1), sym(1/2), sym(1/3), sym(1/4), sym(1/5); sym(1/2), sym(1/3), sym(1/4), sym(1/5), sym(1/6); sym(1/3), sym(1/4), sym(1/5), sym(1/6), sym(1/7); sym(1/4), sym(1/5), sym(1/6), sym(1/7), sym(1/8); sym(1/5), sym(1/6), sym(1/7), sym(1/8), sym(1/9)]

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

d = 1/det(H) 
d = 

-t-2t-32t-43t-54t-65t-74t-83t-92t-1082944-((t - 2)*(t - 3)^2*(t - 4)^3*(t - 5)^4*(t - 6)^5*(t - 7)^4*(t - 8)^3*(t - 9)^2*(t - 10))/82944

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^25/82944 + (25*t^24)/13824 - (5375*t^23)/41472 + (40825*t^22)/6912 - (15940015*t^21)/82944 + (21896665*t^20)/4608 - (240519875*t^19)/2592 + (1268467075*t^18)/864 - (sym('1588946776255')*t^17)/82944 + (sym('2885896606895')*t^16)/13824 - (sym('79493630114675')*t^15)/41472 + (sym('34372691161375')*t^14)/2304 - (sym('8194259295156385')*t^13)/82944 + (sym('7707965729450845')*t^12)/13824 - (sym('55608098247105175')*t^11)/20736 + (sym('37909434298793825')*t^10)/3456 - (sym('197019820623693025')*t^9)/5184 + (sym('10640296363350955')*t^8)/96 - (sym('38821472549340925')*t^7)/144 + (sym('12958201048605475')*t^6)/24 - (sym('1748754621252377')*t^5)/2 + sym('1115685328012530')*t^4 - sym('1078920141906600')*t^3 + sym('742618453752000')*t^2 - sym('323874210240000')*t + sym('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 - 2)*(t - 3)*(t - 4)*(t - 5)*(t - 6)*(t^4 - 18*t^3 + 119*t^2 - 342*t + 360))/576, ((t - 3)*(t - 4)*(t - 5)*(t - 6)*(t - 7)*(t^4 - 17*t^3 + 104*t^2 - 268*t + 240))/144, -((t - 4)*(t - 5)*(t - 6)*(t - 7)*(t - 8)*(t^4 - 16*t^3 + 91*t^2 - 216*t + 180))/96, ((t - 5)*(t - 6)*(t - 7)*(t - 8)*(t - 9)*(t^4 - 15*t^3 + 80*t^2 - 180*t + 144))/144, -((t - 6)*(t - 7)*(t - 8)*(t - 9)*(t - 10)*(t^4 - 14*t^3 + 71*t^2 - 154*t + 120))/576; ((t - 2)*(t - 3)*(t - 4)*(t - 5)*(t - 6)*(t^4 - 22*t^3 + 179*t^2 - 638*t + 840))/144, -((t - 3)*(t - 4)*(t - 5)*(t - 6)*(t - 7)*(t^4 - 21*t^3 + 161*t^2 - 531*t + 630))/36, ((t - 4)*(t - 5)*(t - 6)*(t - 7)*(t - 8)*(t^4 - 20*t^3 + 145*t^2 - 450*t + 504))/24, -((t - 5)*(t - 6)*(t - 7)*(t - 8)*(t - 9)*(t^4 - 19*t^3 + 131*t^2 - 389*t + 420))/36, ((t - 6)*(t - 7)*(t - 8)*(t - 9)*(t - 10)*(t^4 - 18*t^3 + 119*t^2 - 342*t + 360))/144; -((t - 2)*(t - 3)*(t - 4)*(t - 5)*(t - 6)*(t^4 - 26*t^3 + 251*t^2 - 1066*t + 1680))/96, ((t - 3)*(t - 4)*(t - 5)*(t - 6)*(t - 7)*(t^4 - 25*t^3 + 230*t^2 - 920*t + 1344))/24, -((t - 4)*(t - 5)*(t - 6)*(t - 7)*(t - 8)*(t^4 - 24*t^3 + 211*t^2 - 804*t + 1120))/16, ((t - 5)*(t - 6)*(t - 7)*(t - 8)*(t - 9)*(t^4 - 23*t^3 + 194*t^2 - 712*t + 960))/24, -((t - 6)*(t - 7)*(t - 8)*(t - 9)*(t - 10)*(t^4 - 22*t^3 + 179*t^2 - 638*t + 840))/96; ((t - 2)*(t - 3)*(t - 4)*(t - 5)*(t - 6)*(t^4 - 30*t^3 + 335*t^2 - 1650*t + 3024))/144, -((t - 3)*(t - 4)*(t - 5)*(t - 6)*(t - 7)*(t^4 - 29*t^3 + 311*t^2 - 1459*t + 2520))/36, ((t - 4)*(t - 5)*(t - 6)*(t - 7)*(t - 8)*(t^4 - 28*t^3 + 289*t^2 - 1302*t + 2160))/24, -((t - 5)*(t - 6)*(t - 7)*(t - 8)*(t - 9)*(t^4 - 27*t^3 + 269*t^2 - 1173*t + 1890))/36, ((t - 6)*(t - 7)*(t - 8)*(t - 9)*(t - 10)*(t^4 - 26*t^3 + 251*t^2 - 1066*t + 1680))/144; -((t - 2)*(t - 3)*(t - 4)*(t - 5)*(t - 6)*(t^4 - 34*t^3 + 431*t^2 - 2414*t + 5040))/576, ((t - 3)*(t - 4)*(t - 5)*(t - 6)*(t - 7)*(t^4 - 33*t^3 + 404*t^2 - 2172*t + 4320))/144, -((t - 4)*(t - 5)*(t - 6)*(t - 7)*(t - 8)*(t^4 - 32*t^3 + 379*t^2 - 1968*t + 3780))/96, ((t - 5)*(t - 6)*(t - 7)*(t - 8)*(t - 9)*(t^4 - 31*t^3 + 356*t^2 - 1796*t + 3360))/144, -((t - 6)*(t - 7)*(t - 8)*(t - 9)*(t - 10)*(t^4 - 30*t^3 + 335*t^2 - 1650*t + 3024))/576]

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

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

(25-3001050-1400630-3004800-1890026880-126001050-1890079380-11760056700-140026880-117600179200-88200630-1260056700-8820044100)[sym(25), -sym(300), sym(1050), -sym(1400), sym(630); -sym(300), sym(4800), -sym(18900), sym(26880), -sym(12600); sym(1050), -sym(18900), sym(79380), -sym(117600), sym(56700); -sym(1400), sym(26880), -sym(117600), sym(179200), -sym(88200); sym(630), -sym(12600), sym(56700), -sym(88200), sym(44100)]

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)[-sym(9), sym(11), -sym(21), sym(63), -sym(252); sym(70), -sym(69), sym(141), -sym(421), sym(1684); -sym(575), sym(575), -sym(1149), sym(3451), -sym(13801); sym(3891), -sym(3891), sym(7782), -sym(23345), sym(93365); sym(1024), -sym(1024), sym(2048), -sym(6144), sym(24572)]

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

A^5 
ans = 

(0000000000000000000000000)[sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), sym(0)]

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

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

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

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

lambda = eig(A) 
lambda = 

(00000)[sym(0); sym(0); sym(0); sym(0); sym(0)]

数値計算は丸め誤差を伴い、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)[vpa('0.0005617448486395847'); vpa('0.0001737348850386136') - vpa('0.0005342985684139864i'); vpa('0.0001737348850386136') + vpa('0.0005342985684139864i'); - vpa('0.000454607309358406') + vpa('0.0003303865815979566i'); - vpa('0.000454607309358406') - vpa('0.0003303865815979566i')]

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

J = jordan(A) 
J = 

(0100000100000100000100000)[sym(0), sym(1), sym(0), sym(0), sym(0); sym(0), sym(0), sym(1), sym(0), sym(0); sym(0), sym(0), sym(0), sym(1), sym(0); sym(0), sym(0), sym(0), sym(0), sym(1); sym(0), sym(0), sym(0), sym(0), sym(0)]

行列の指数関数 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)[- (2*t^3)/3 + (11*t^2)/2 - 9*t + 1, (t*(4*t^2 - 27*t + 33))/3, -(t*(20*t^2 - 117*t + 126))/6, (t*(32*t^2 - 174*t + 189))/3, -(t*(85*t^2 - 464*t + 504))/2; -(t*(7*t^3 - 81*t^2 + 230*t - 140))/2, 7*t^4 - 67*t^3 + (301*t^2)/2 - 69*t + 1, -(t*(35*t^3 - 293*t^2 + 598*t - 282))/2, (t*(112*t^3 - 876*t^2 + 1799*t - 842))/2, -(t*(1785*t^3 - 14012*t^2 + 28776*t - 13472))/8; (t*(142*t^3 - 1710*t^2 + 5151*t - 3450))/6, -(t*(142*t^3 - 1426*t^2 + 3438*t - 1725))/3, (355*t^4)/3 - (3139*t^3)/3 + (4585*t^2)/2 - 1149*t + 1, -(t*(1136*t^3 - 9420*t^2 + 20625*t - 10353))/3, (t*(18105*t^3 - 150646*t^2 + 329952*t - 165612))/12; -(t*(973*t^3 - 11675*t^2 + 35022*t - 23346))/6, (t*(1946*t^3 - 19458*t^2 + 46695*t - 23346))/6, -(t*(4865*t^3 - 42807*t^2 + 93390*t - 46692))/6, (7784*t^4)/3 - (64210*t^3)/3 + (93391*t^2)/2 - 23345*t + 1, -(t*(248115*t^3 - 2053748*t^2 + 4482036*t - 2240760))/24; -(128*t*(t^3 - 12*t^2 + 36*t - 24))/3, (256*t*(t^3 - 10*t^2 + 24*t - 12))/3, -(128*t*(5*t^3 - 44*t^2 + 96*t - 48))/3, (512*t*(4*t^3 - 33*t^2 + 72*t - 36))/3, - 2720*t^4 + (67552*t^3)/3 - 49144*t^2 + 24572*t + 1]

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

X = exp(t*A) 
X = 

(e-9te11te-21te63te-252te70te-69te141te-421te1684te-575te575te-1149te3451te-13801te3891te-3891te7782te-23345te93365te1024te-1024te2048te-6144te24572t)[exp(-9*t), exp(11*t), exp(-21*t), exp(63*t), exp(-252*t); exp(70*t), exp(-69*t), exp(141*t), exp(-421*t), exp(1684*t); exp(-575*t), exp(575*t), exp(-1149*t), exp(3451*t), exp(-13801*t); exp(3891*t), exp(-3891*t), exp(7782*t), exp(-23345*t), exp(93365*t); exp(1024*t), exp(-1024*t), exp(2048*t), exp(-6144*t), exp(24572*t)]