Символьные матричные расчеты

Этот пример показывает, как выполнить простые матричные расчеты с помощью 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.567050691098230795533011005520724633949315252334'); vpa ('0.2085342186110133359050025100688200550385820260343'); vpa ('0.011407491623419806559451458886658934504234843052664'); vpa ('0.00030589804015119172687994978406922282565614514909247'); 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)]

Эта матрица является «nilpotent». Это пятая степень - нулевая матрица.

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 = λ5лямбда ^ sym (5)

Теперь вы сможете вычислить собственные значения матрицы в голове. Они являются нулями уравнения lambda ^ 5 = 0.

Символьные расчеты могут точно найти собственные значения.

lambda = eig(A) 
lambda = 

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

Численное вычисление включает округлую ошибку и находит нули уравнения, которое является чем-то вроде лямбды ^ 5 = eps * norm (A) Таким образом, вычисленные собственные значения примерно лямбда = (eps * norm (A)) ^ (1/5) Вот собственные значения, вычисленные Symbolic Toolbox с помощью Не очевидно, что все они должны быть нулевыми.

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 (лямбда (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 - 1500,)/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)]