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/315*x^4+735781/2116800*x^3-852401/222264000*x^2+61501/53343360000*x-1/266716800000
Try to factor the characteristic polynomial.
factor(p)
ans = x^5-563/315*x^4+735781/2116800*x^3-852401/222264000*x^2+61501/53343360000*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 =
.32879287721718629571150047605447313997367890616371e-5
.30589804015119172687949784069272282565614514909017e-3
.11407491623419806559451458866589345042348430526627e-1
.20853421861101333590500251006882005503858202260342
1.5670506910982307955330110055207246339493152522336
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/(2-t), 1/(3-t), 1/(4-t), 1/(5-t), 1/(6-t)] [ 1/(3-t), 1/(4-t), 1/(5-t), 1/(6-t), 1/(7-t)] [ 1/(4-t), 1/(5-t), 1/(6-t), 1/(7-t), 1/(8-t)] [ 1/(5-t), 1/(6-t), 1/(7-t), 1/(8-t), 1/(9-t)] [ 1/(6-t), 1/(7-t), 1/(8-t), 1/(9-t), 1/(10-t)]
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 =
-1/82944*(-2+t)*(-4+t)^3*(-6+t)^5*(-8+t)^3*(-10+t)*(-9+t)^2*(-7+t)^4*(-5+t)^4*(-3+t)^2
d =
-323874210240000*t+742618453752000*t^2-1078920141906600*t^3+1115685328012530*t^4-1748754621252377/2*t^5+12958201048605475/24*t^6-38821472549340925/144*t^7+10640296363350955/96*t^8-197019820623693025/5184*t^9+37909434298793825/3456*t^10-55608098247105175/20736*t^11+7707965729450845/13824*t^12-8194259295156385/82944*t^13+34372691161375/2304*t^14-79493630114675/41472*t^15+2885896606895/13824*t^16-1588946776255/82944*t^17+1268467075/864*t^18-240519875/2592*t^19+21896665/4608*t^20-15940015/82944*t^21+40825/6912*t^22-5375/41472*t^23+67212633600000+25/13824*t^24-1/82944*t^25
2 3
-323874210240000 t + 742618453752000 t - 1078920141906600 t
4 5 12958201048605475 6
+ 1115685328012530 t - 1748754621252377/2 t + ----------------- t
24
38821472549340925 7 10640296363350955 8 197019820623693025 9
- ----------------- t + ----------------- t - ------------------ t
144 96 5184
37909434298793825 10 55608098247105175 11
+ ----------------- t - ----------------- t
3456 20736
7707965729450845 12 8194259295156385 13 34372691161375 14
+ ---------------- t - ---------------- t + -------------- t
13824 82944 2304
79493630114675 15 2885896606895 16 1588946776255 17
- -------------- t + ------------- t - ------------- t
41472 13824 82944
1268467075 18 240519875 19 21896665 20 15940015 21
+ ---------- t - --------- t + -------- t - -------- t
864 2592 4608 82944
40825 22 5375 23 25 24 25
+ ----- t - ----- t + 67212633600000 + ----- t - 1/82944 t
6912 41472 13824
The elements of the inverse are also polynomials in t.
X = inv(H)
X = [ -1/576*(-4+t)^2*(-6+t)^2*(-5+t)^2*(-2+t)*(-3+t)^2, 1/144*(-3+t)*(-6+t)^2*(-7+t)*(-5+t)^2*(-4+t)^2*(-2+t), -1/96*(-3+t)*(-5+t)^2*(-8+t)*(-7+t)*(-6+t)^2*(-4+t)*(-2+t), 1/144*(-3+t)*(-5+t)*(-7+t)*(-8+t)*(-9+t)*(-6+t)^2*(-4+t)*(-2+t), -1/576*(-3+t)*(-5+t)*(-7+t)*(-9+t)*(-8+t)*(-6+t)*(-4+t)*(-2+t)*(-10+t)] [ 1/144*(-3+t)*(-6+t)^2*(-7+t)*(-5+t)^2*(-4+t)^2*(-2+t), -1/36*(-6+t)^2*(-7+t)^2*(-4+t)*(-5+t)^2*(-3+t)^2, 1/24*(-5+t)*(-8+t)*(-7+t)^2*(-6+t)^2*(-4+t)^2*(-3+t), -1/36*(-3+t)*(-6+t)*(-9+t)*(-7+t)^2*(-8+t)*(-5+t)^2*(-4+t), 1/144*(-5+t)*(-7+t)*(-9+t)*(-8+t)*(-6+t)^2*(-4+t)*(-3+t)*(-10+t)] [ -1/96*(-3+t)*(-5+t)^2*(-8+t)*(-7+t)*(-6+t)^2*(-4+t)*(-2+t), 1/24*(-5+t)*(-8+t)*(-7+t)^2*(-6+t)^2*(-4+t)^2*(-3+t), -1/16*(-4+t)^2*(-8+t)^2*(-6+t)*(-7+t)^2*(-5+t)^2, 1/24*(-4+t)*(-7+t)*(-8+t)^2*(-9+t)*(-5+t)^2*(-6+t)^2, -1/96*(-4+t)*(-7+t)^2*(-9+t)*(-8+t)*(-6+t)^2*(-5+t)*(-10+t)] [ 1/144*(-3+t)*(-5+t)*(-7+t)*(-8+t)*(-9+t)*(-6+t)^2*(-4+t)*(-2+t), -1/36*(-3+t)*(-6+t)*(-9+t)*(-7+t)^2*(-8+t)*(-5+t)^2*(-4+t), 1/24*(-4+t)*(-7+t)*(-8+t)^2*(-9+t)*(-5+t)^2*(-6+t)^2, -1/36*(-6+t)^2*(-8+t)*(-5+t)^2*(-7+t)^2*(-9+t)^2, 1/144*(-5+t)*(-7+t)^2*(-10+t)*(-8+t)^2*(-9+t)*(-6+t)^2] [ -1/576*(-3+t)*(-5+t)*(-7+t)*(-9+t)*(-8+t)*(-6+t)*(-4+t)*(-2+t)*(-10+t), 1/144*(-5+t)*(-7+t)*(-9+t)*(-8+t)*(-6+t)^2*(-4+t)*(-3+t)*(-10+t), -1/96*(-4+t)*(-7+t)^2*(-9+t)*(-8+t)*(-6+t)^2*(-5+t)*(-10+t), 1/144*(-5+t)*(-7+t)^2*(-10+t)*(-8+t)^2*(-9+t)*(-6+t)^2, -1/576*(-8+t)^2*(-7+t)^2*(-6+t)^2*(-10+t)*(-9+t)^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 =
-.2714743126973639e-2
-.841289961972077e-3+.2582690715993584e-2*i
-.841289961972077e-3-.2582690715993584e-2*i
.219866152544703e-2+.1598907770696703e-2*i
.219866152544703e-2-.1598907770696703e-2*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 = expm(t*A)
E = [ 1-9*t+11/2*t^2-2/3*t^3, 4/3*t^3-9*t^2+11*t, -10/3*t^3+39/2*t^2-21*t, 32/3*t^3-58*t^2+63*t, -85/2*t^3+232*t^2-252*t] [ 70*t-115*t^2+81/2*t^3-7/2*t^4, 1+7*t^4-67*t^3+301/2*t^2-69*t, -35/2*t^4+293/2*t^3-299*t^2+141*t, 56*t^4-438*t^3+1799/2*t^2-421*t, -1785/8*t^4+3503/2*t^3-3597*t^2+1684*t] [ -575*t+1717/2*t^2-285*t^3+71/3*t^4, -142/3*t^4+1426/3*t^3-1146*t^2+575*t, 1+355/3*t^4-3139/3*t^3+4585/2*t^2-1149*t, -1136/3*t^4+3140*t^3-6875*t^2+3451*t, 6035/4*t^4-75323/6*t^3+27496*t^2-13801*t] [ 3891*t-5837*t^2+11675/6*t^3-973/6*t^4, 973/3*t^4-3243*t^3+15565/2*t^2-3891*t, -4865/6*t^4+14269/2*t^3-15565*t^2+7782*t, 1+7784/3*t^4-64210/3*t^3+93391/2*t^2-23345*t, -82705/8*t^4+513437/6*t^3-373503/2*t^2+93365*t] [ 1024*t-1536*t^2+512*t^3-128/3*t^4, 256/3*t^4-2560/3*t^3+2048*t^2-1024*t, -640/3*t^4+5632/3*t^3-4096*t^2+2048*t, 2048/3*t^4-5632*t^3+12288*t^2-6144*t, 1-2720*t^4+67552/3*t^3-49144*t^2+24572*t]
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)]