Linear Algebra and Its Applications, Exercise 1.7.8

Exercise 1.7.8. Take the 10 by 10 Hilbert matrix A for which a_{ij} = 1/(i+j-1) and solve the equation Ax = (1, 0, \dotsc, 0). Make a minor change to A or b and see how the solution x changes.

Answer: We can use the open source statistics software R for this exercise. The following R commands will define a 10 by 10 Hilbert matrix. The first command defines A as a 10 by 10 matrix with all values set to 0, and the second command sets each entry a_{ij} to the proper value for a Hilbert matrix. (Note that there may be a faster or more elegant way to do this in R, but this way is easily understandable.)

> A <- array(0, dim=c(10,10))
> for (i in 1:10) {for (j in 1:10) {A[i, j] <- 1/(i+j-1)}}
>

We then print the matrix to double-check that it appears to be correct (we use the option() function to limit the number of digits printed per entry and increase the width of the printable area, so that the entire matrix can be displayed more compactly):

> options(digits=3)
> options(width=90)
> A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1.000 0.500 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000
[2,] 0.500 0.333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909
[3,] 0.333 0.250 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833
[4,] 0.250 0.200 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769
[5,] 0.200 0.167 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714
[6,] 0.167 0.143 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667
[7,] 0.143 0.125 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625
[8,] 0.125 0.111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588
[9,] 0.111 0.100 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556
[10,] 0.100 0.091 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526
>

We next define b = (1, 0, \dotsc, 0):

> b <- array(0, dim=c(10))
> b[1] <- 1
> b
 [1] 1 0 0 0 0 0 0 0 0 0
>

We can now use the R solve() function to solve for x:

> x <- solve(A, b)
> x
[1]      100    -4950    79195  -600553  2522294 -6305679  9608580 -8750614  4375282
[10]  -923666
>

We then multiply A by x to verify that x is indeed a solution to Ax = b:

> A %*% x
           [,1]
 [1,]  1.00e+00
 [2,] -1.89e-10
 [3,] -2.04e-10
 [4,] -1.60e-10
 [5,] -1.02e-10
 [6,]  3.64e-11
 [7,]  2.18e-11
 [8,]  7.28e-11
 [9,]  0.00e+00
[10,]  0.00e+00
> 

Finally we tweak b very slightly and see how the solution to Ax = b changes:

> b2 <- b
> b2[10] <- .01
> b2
[1] 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
> x2 <- solve(A, b2)
> x2
[1] -9.14e+03  8.26e+05 -1.82e+07  1.70e+08 -8.30e+08  2.32e+09 -3.87e+09  3.80e+09
[9] -2.02e+09  4.48e+08
>

Note that the very minor change to b caused the entries of the solution x to change by one or more orders of magnitude.

NOTE: This continues a series of posts containing worked out exercises from the (out of print) book Linear Algebra and Its Applications, Third Edition by Gilbert Strang.

If you find these posts useful I encourage you to also check out the more current Linear Algebra and Its Applications, Fourth Edition, Dr Strang’s introductory textbook Introduction to Linear Algebra, Fourth Edition and the accompanying free online course, and Dr Strang’s other books.

This entry was posted in linear algebra. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s