Store 
Symbolic Math Toolbox 5.3
Symbolic Matrix Computation
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.
p = poly(H)
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 + 30 7505*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.567050691098230795533011005520724633949315252233401
0.2085342186110133359050025100688200550385820226034302
0.01140749162341980655945145886658934504234843052664014
0.0003058980401511917268794978406927228256561451490924693
0.000003287928772171862957115004760544731399736789023074625
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.0000 0.5000 0.3333 0.2500 0.2000
0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.2500 0.2000 0.1667 0.1429
0.2500 0.2000 0.1667 0.1429 0.1250
0.2000 0.1667 0.1429 0.1250 0.1111
The reciprocal of the determinant is a polynomial in t.
d = 1/det(H) d = expand(d) pretty(d)
d =
- t^25/82944 + (25*t^24)/13824 - (5375*t^23)/41472 + (40825*t^22)/6912 - (15
940015*t^21)/82944 + (21896665*t^20)/4608 - (240519875*t^19)/2592 + (1268467
075*t^18)/864 - (1588946776255*t^17)/82944 + (2885896606895*t^16)/13824 - (7
9493630114675*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 + (106402963
63350955*t^8)/96 - (38821472549340925*t^7)/144 + (12958201048605475*t^6)/24
- (1748754621252377*t^5)/2 + 1115685328012530*t^4 - 1078920141906600*t^3 + 7
42618453752000*t^2 - 323874210240000*t + 67212633600000
d =
- t^25/82944 + (25*t^24)/13824 - (5375*t^23)/41472 + (40825*t^22)/6912 - (15
940015*t^21)/82944 + (21896665*t^20)/4608 - (240519875*t^19)/2592 + (1268467
075*t^18)/864 - (1588946776255*t^17)/82944 + (2885896606895*t^16)/13824 - (7
9493630114675*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 + (106402963
63350955*t^8)/96 - (38821472549340925*t^7)/144 + (12958201048605475*t^6)/24
- (1748754621252377*t^5)/2 + 1115685328012530*t^4 - 1078920141906600*t^3 + 7
42618453752000*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/576 - 1/288)*(t^4 - 18*t^3 + 119*t^2 - 342*t + 360)^2,
(t^3 - 15*t^2 + 74*t - 120)^2
*(t^3/144 - t^2/12 + (41*t)/144 - 7/24), -(t^2 - 11*t + 30)^2*(t^5/96
- t^4/4 + (217*t^3)/96 - (153*t^2)/16 + (227*t)/12 - 14), (t - 6)^2*(t^7/144
- (19*t^6)/72 + (299*t^5)/72 - (1259*t^4)/36 + (24433*t^3)/144 - (34039*t^2
)/72 + (4189*t)/6 - 420), - t^9/576 + (3*t^8)/32 - (211*t^7)/96 + (469*t^6)/
16 - (46963*t^5)/192 + (42287*t^4)/32 - (663941*t^3)/144 + (79913*t^2)/8 - (
24305*t)/2 + 6300]
[
(t^3 - 15*t^2 + 74*t - 120)^2*(t^3/144 - t^2/12 + (41*t)/144 - 7/24),
-(t/36 - 1/9)*(
t^4 - 21*t^3 + 161*t^2 - 531*t + 630)^2, (t^3 - 17*
t^2 + 94*t - 168)^2*(t^3/24 - (2*t^2)/3 + (79*t)/24 - 5),
-(t^2 - 12*t + 35)^2*(t^5/36 - (5*t^4)/6 + (347*t^3)/36 - (1
07*t^2)/2 + 142*t - 144), (t - 6)^2*(t^7/144 - (23*t^
6)/72 + (443*t^5)/72 - (2311*t^4)/36 + (56305*t^3)/144 - (99899*t^2)/72 + (1
5905*t)/6 - 2100)]
[ -(t^2 - 11*t + 3
0)^2*(t^5/96 - t^4/4 + (217*t^3)/96 - (153*t^2)/16 + (227*t)/12 - 14),
(t^3 - 17*t^2 + 94*t - 168)
^2*(t^3/24 - (2*t^2)/3 + (79*t)/24 - 5),
-(t/16 - 3/8)*(t^4 - 24*t^3 + 211*t^2 - 804*t + 1120)^2,
(t^3 - 19*t^2 + 118*t - 240)^2*(t^3/24 - (5*t^2
)/6 + (127*t)/24 - 21/2),
-(t^2 - 13*t + 42)^2*(t^5/96 - (3*t^4)/8 + (505*t^3)/96 - (573*t^2)/16 + (
1415*t)/12 - 150)]
[ (t - 6)^2*(t^7/144 - (19*t^6)/72 + (299*t^5)/72 -
(1259*t^4)/36 + (24433*t^3)/144 - (34039*t^2)/72 + (4189*t)/6 - 420),
-(t^2 - 12*t + 35)^2*(t^5/36 - (5*t^4)/6 + (3
47*t^3)/36 - (107*t^2)/2 + 142*t - 144), (t^3 - 19*t^2 +
118*t - 240)^2*(t^3/24 - (5*t^2)/6 + (127*t)/24 - 21/2),
-(t/36 - 2/9)*(t^4 - 27*t^3 + 26
9*t^2 - 1173*t + 1890)^2,
(t^3 - 21*t^2 + 146*t - 336)^2*(t^3/144 - t^2/6 + (1
85*t)/144 - 25/8)]
[ - t^9/576 + (3*t^8)/32 - (211*t^7)/96 + (469*t^6)/16 - (46963*t^5)/192 + (
42287*t^4)/32 - (663941*t^3)/144 + (79913*t^2)/8 - (24305*t)/2 + 6300, (t -
6)^2*(t^7/144 - (23*t^6)/72 + (443*t^5)/72 - (2311*t^4)/36 + (56305*t^3)/144
- (99899*t^2)/72 + (15905*t)/6 - 2100), -(t^2 - 13*t + 42)^2*(t^5/96 - (3*t
^4)/8 + (505*t^3)/96 - (573*t^2)/16 + (1415*t)/12 - 150),
(t^3 - 21*t^2 + 146*t - 336)^2*(t^3/144 - t^2
/6 + (185*t)/144 - 25/8),
-(t/576 - 5/288)*(t^4 - 30*t^3 + 335*t^2 -
1650*t + 3024)^2]
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 = poly(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.000561744847392736596
0.000534298567227930807*i + 0.00017373488465260121
0.00017373488465260121 - 0.0005342985672279308*i
0.000330386580864093695*i - 0.0004546073083489695
- 0.0003303865808640937*i - 0.0004546073083489695
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*(2*t - 3)*(16*t - 63))/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 + (3
01*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 = [ 1/exp(9*t), exp(11*t), 1/exp(21*t), exp(63*t), 1/exp(252*t) ] [ exp(70*t), 1/exp(69*t), exp(141*t), 1/exp(421*t), exp(1684*t) ] [ 1/exp(575*t), exp(575*t), 1/exp(1149*t), exp(3451*t), 1/exp(13801*t) ] [ exp(3891*t), 1/exp(3891*t), exp(7782*t), 1/exp(23345*t), exp(93365*t) ] [ exp(1024*t), 1/exp(1024*t), exp(2048*t), 1/exp(6144*t), exp(24572*t) ]
