Расчет символьной матрицы

В этом примере показано, как выполнить простые матричные расчеты с помощью 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') т + 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)]

Эта матрица является "нильпотентной". Это - пятая степень, нулевая матрица.

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), Вот собственные значения, вычисленные Символьным Тулбоксом с помощью 16 арифметик плавающей точки цифры. Не очевидно, что они должны все быть нулем.

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)]