Homework 1

We consider a system of linear equations A x = b , where A , b are given matrices nxn and nx1.

Contents

Problem 1

Find x for

A = [1 1 1 1 1 1;1 2 3 4 5 6;1 3 6 10 15 21;1 4 10 20 35 56;1 5 15 35 70 126;1 6 21 56 126 252]
b = [1 2 3 4 5 6]'
A =

     1     1     1     1     1     1
     1     2     3     4     5     6
     1     3     6    10    15    21
     1     4    10    20    35    56
     1     5    15    35    70   126
     1     6    21    56   126   252


b =

     1
     2
     3
     4
     5
     6

Problem 2

Find x for

A = [1 2 3;4 5 6;7 8 9]
b = [6 15 20]'
A =

     1     2     3
     4     5     6
     7     8     9


b =

     6
    15
    20

Problem 3

Find x for

A = [1 2 3 4;4 5 6 7;5 7 9 11;3 3 3 3]
b = [30 60 90 30]'
A =

     1     2     3     4
     4     5     6     7
     5     7     9    11
     3     3     3     3


b =

    30
    60
    90
    30

Solving systems of linear equations in Matlab

In the lecture, we all know that if the matrix A is nonsingular, we have a unique solution, but if A is singular, we have either no solution or infinitely many solutions. Therefore, the fist thing we want to know is wherther A is singular or nonsingular.

Checking singularity

In all textbooks, we are told that if det A = 0 then A is singular. However, in Matlab, using determiants to check singularity is not reliable because of floating-point arithmetic operations https://en.wikipedia.org/wiki/Floating_point

For example

A = [1 2 3;4 5 6;7 8 9]
detA = det(A) %calculate determinant of A
A =

     1     2     3
     4     5     6
     7     8     9


detA =

   6.6613e-16

Matrix A is singular, but det A isn't zero. If you are using Matlab 2010a or older version, the result may be zero but for newer version (I'm using Malab 2014b), the result is not zero.

The realible way to check singularity of a square nxn matrix A is to calculate rank of A. If rank A < n, the matrix is singular.

A = [1 2 3;4 5 6;7 8 9]
rankA = rank(A) %calculate rank of A
A =

     1     2     3
     4     5     6
     7     8     9


rankA =

     2

A is a 3x3 matrix, rank A = 2 < 3, so A is singular

Nonsingular matrix

If A is nonsingular, the equation A*x* = b has unique solution x = inv(A)*b or A\b

For example

A = [8 1 6;3 5 7;4 9 2]
b = [1 2 3]'
x = inv(A)*b
% or
x = A\b
A =

     8     1     6
     3     5     7
     4     9     2


b =

     1
     2
     3


x =

    0.0500
    0.3000
    0.0500


x =

    0.0500
    0.3000
    0.0500

Singular matrix

If A is singular, the solution to Ax = b either does not exist, or is not unique. You can determine whether Ax = b has a solution by finding the row reduced echelon form of the augmented matrix [A b] https://en.wikipedia.org/wiki/Row_echelon_form . To do so for this example, enter rref([A b]). If there is a row which contains all zeros except for the last entry, the equation does not have a solution.

For example

A = [1 1;1 1]
b = [1 2]'
rref([A b])
A =

     1     1
     1     1


b =

     1
     2


ans =

     1     1     0
     0     0     1

The bottom row has all zero elements except one entry therefore the equation Ax = b does not have a solution.

If the equation Ax = b has solutions. The general solution to a system of linear equations Ax = b describes all possible solutions. You can find the general solution by:

1. Solving the corresponding homogeneous system Ax = 0. Do this using the null command, by typing null(A). This returns a basis for the solution space to Ax = 0. Any solution is a linear combination of basis vectors.

2. Finding a particular solution to the nonhomogeneous system Ax = b by using the command pinv(A)*b.

Example

A = [1 2 3;2 4 6;3 6 9]
b = [1 2 3]'
rref([A b])
A =

     1     2     3
     2     4     6
     3     6     9


b =

     1
     2
     3


ans =

     1     2     3     1
     0     0     0     0
     0     0     0     0

From the result we can conclude that the equation has solutions.

Find vectors span the null space of A

null(A)
ans =

    0.6777    0.6850
    0.4872   -0.6906
   -0.5507    0.2320

Because of all elements are rational elements, we can find the vectors in rational numbers

null(A,'r') % You can use sym(null(A))
ans =

    -2    -3
     1     0
     0     1

Any linear combination of two vectors (-2 1 0)' and (-3 0 1) is a soltion of Ax = 0;

Next, we find a particular solution of Ax = b by a command pinv(A)*b

x_p = pinv(A)*b
x_p =

    0.0714
    0.1429
    0.2143

we can convert x_p to rational numbers

sym(x_p)
 
ans =
 
 1/14
  1/7
 3/14
 

The general soltion of the equation Ax = b is x1*[-2 1 0]' + x2*[-3 0 1]' + [1/14 1/7 3/14]', where x1, x2 are two arbitrary numbers.