Symbolic Matrix Computation

This example shows how to perform simple matrix computations using Symbolic Math Toolbox™.

Generate a possibly familiar test matrix, the 5-by-5 Hilbert matrix.

H = sym(hilb(5)) 
H = 

(1121314151213141516131415161714151617181516171819)

The determinant is very small.

d = det(H) 
d = 

1266716800000

The elements of the inverse are integers.

X = inv(H) 
X = 

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

Verify that the inverse is correct.

I = X*H
I = 

(1000001000001000001000001)

Find the characteristic polynomial.

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

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

Try to factor the characteristic polynomial.

factor(p) 
ans = 

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

The result indicates that the characteristic polynomial cannot be factored over the rational numbers.

Compute the 50 digit numerical approximations to the eigenvalues.

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

(1.56705069109823079553301100552072463394931525223340.208534218611013335905002510068820055038582022603430.011407491623419806559451458866589345042348430526640.000305898040151191726879497840692722825656145149092470.0000032879287721718629571150047605447313997367890230746)

Create a generalized Hilbert matrix involving a free variable, 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

Substituting t=1 retrieves the original Hilbert matrix.

subs(H,t,1) 
ans = 

(1121314151213141516131415161714151617181516171819)

The reciprocal of the determinant is a polynomial in 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

The elements of the inverse are also polynomials in 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

Substituting t=1 generates the Hilbert inverse.

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

Investigate a different example.

A = sym(gallery(5)) 
A = 

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

This matrix is "nilpotent". It's fifth power is the zero matrix.

A^5 
ans = 

(0000000000000000000000000)

Because this matrix is nilpotent, its characteristic polynomial is very simple.

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

You should now be able to compute the matrix eigenvalues in your head. They are the zeros of the equation lambda^5 = 0.

Symbolic computation can find the eigenvalues exactly.

lambda = eig(A) 
lambda = 

(00000)

Numeric computation involves roundoff error and finds the zeros of an equation that is something like lambda^5 = eps*norm(A) So the computed eigenvalues are roughly lambda = (eps*norm(A))^(1/5) Here are the eigenvalues, computed by the Symbolic Toolbox using 16 digit floating point arithmetic. It is not obvious that they should all be zero.

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

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

This matrix is also "defective". It is not similar to a diagonal matrix. Its Jordan Canonical Form is not diagonal.

J = jordan(A) 
J = 

(0100000100000100000100000)

The matrix exponential, expm(t*A), is usually expressed in terms of scalar exponentials involving the eigenvalues, exp(lambda(i)*t). But for this matrix, the elements of expm(t*A) are all polynomials in 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)

By the way, the function "exp" computes element-by-element exponentials.

X = exp(t*A) 
X = 

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