Quantcast

Documentation Center

  • Trial Software
  • Product Updates

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 =
 
[   1, 1/2, 1/3, 1/4, 1/5]
[ 1/2, 1/3, 1/4, 1/5, 1/6]
[ 1/3, 1/4, 1/5, 1/6, 1/7]
[ 1/4, 1/5, 1/6, 1/7, 1/8]
[ 1/5, 1/6, 1/7, 1/8, 1/9]
 

The determinant is very small.

d = det(H)
 
d =
 
1/266716800000
 

The elements of the inverse are integers.

X = inv(H)
 
X =
 
[    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]
 

Verify that the inverse is correct.

I = X*H
 
I =
 
[ 1, 0, 0, 0, 0]
[ 0, 1, 0, 0, 0]
[ 0, 0, 1, 0, 0]
[ 0, 0, 0, 1, 0]
[ 0, 0, 0, 0, 1]
 

Find the characteristic polynomial.

syms x;
p = charpoly(H,x)
 
p =
 
x^5 - (563*x^4)/315 + (735781*x^3)/2116800 - (852401*x^2)/222264000 + (61501*x)/53343360000 - 1/266716800000
 

Try to factor the characteristic polynomial.

factor(p)
 
ans =
 
(266716800000*x^5 - 476703360000*x^4 + 92708406000*x^3 - 1022881200*x^2 + 307505*x - 1)/266716800000
 

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

Compute 50 digit numerical approximations to the eigenvalues.

digits(50)
e = eig(vpa(H))
 
e =
 
       1.5670506910982307955330110055207246339493152522334
      0.20853421861101333590500251006882005503858202260343
      0.01140749162341980655945145886658934504234843052664
   0.00030589804015119172687949784069272282565614514909247
 0.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 =
 
[ -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)]
 

Substituting t = 1 retrieves the original Hilbert matrix.

subs(H,t,1)
 
ans =
 
[   1, 1/2, 1/3, 1/4, 1/5]
[ 1/2, 1/3, 1/4, 1/5, 1/6]
[ 1/3, 1/4, 1/5, 1/6, 1/7]
[ 1/4, 1/5, 1/6, 1/7, 1/8]
[ 1/5, 1/6, 1/7, 1/8, 1/9]
 

The reciprocal of the determinant is a polynomial in t.

d = 1/det(H)

d = expand(d)

pretty(d)
 
d =
 
-((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 =
 
- 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 - (1588946776255*t^17)/82944 + (2885896606895*t^16)/13824 - (79493630114675*t^15)/41472 + (34372691161375*t^14)/2304 - (8194259295156385*t^13)/82944 + (7707965729450845*t^12)/13824 - (55608098247105175*t^11)/20736 + (37909434298793825*t^10)/3456 - (197019820623693025*t^9)/5184 + (10640296363350955*t^8)/96 - (38821472549340925*t^7)/144 + (12958201048605475*t^6)/24 - (1748754621252377*t^5)/2 + 1115685328012530*t^4 - 1078920141906600*t^3 + 742618453752000*t^2 - 323874210240000*t + 67212633600000
 
    25        24         23          22             21             20
   t      25 t     5375 t     40825 t     15940015 t     21896665 t
- ----- + ------ - -------- + --------- - ------------ + ------------
  82944    13824     41472       6912         82944          4608

                19               18                  17                  16
     240519875 t     1268467075 t     1588946776255 t     2885896606895 t
   - ------------- + -------------- - ----------------- + -----------------
          2592             864              82944               13824

                     15                   14                     13
     79493630114675 t     34372691161375 t     8194259295156385 t
   - ------------------ + ------------------ - --------------------
            41472                2304                  82944

                       12                      11                      10
     7707965729450845 t     55608098247105175 t     37909434298793825 t
   + -------------------- - --------------------- + ---------------------
             13824                  20736                    3456

                         9                      8                      7
     197019820623693025 t    10640296363350955 t    38821472549340925 t
   - --------------------- + -------------------- - --------------------
              5184                    96                     144

                        6                     5
     12958201048605475 t    1748754621252377 t                      4
   + -------------------- - ------------------- + 1115685328012530 t
              24                     2

                       3                    2
   - 1078920141906600 t  + 742618453752000 t  - 323874210240000 t

   + 67212633600000

The elements of the inverse are also polynomials in t.

X = inv(H)
 
X =
 
[   -((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]
 

Substituting t = 1 generates the Hilbert inverse.

X = subs(X,t,'1')
X = double(X)
 
X =
 
[    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]
 

X =

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

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

A^5
 
ans =
 
[ 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0]
 

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

p = charpoly(A,'lambda')
 
p =
 
lambda^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 =
 
 0
 0
 0
 0
 0
 

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.0005617448486395847
  0.0001737348850386136 - 0.0005342985684139864*i
  0.0001737348850386136 + 0.0005342985684139864*i
 - 0.000454607309358406 + 0.0003303865815979566*i
 - 0.000454607309358406 - 0.0003303865815979566*i
 

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

J = jordan(A)
 
J =
 
[ 0, 1, 0, 0, 0]
[ 0, 0, 1, 0, 0]
[ 0, 0, 0, 1, 0]
[ 0, 0, 0, 0, 1]
[ 0, 0, 0, 0, 0]
 

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 =
 
[             - (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]
 

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

X = exp(t*A)
 
X =
 
[   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)]
 
Was this topic helpful?