Solving N Linear Equations for N Unknowns With MATLAB

 

 

How to do it

 

·       In a wide variety of fields it is necessary to solve N equations for N unknowns.  This is required in electric circuit analysis, weather models, and even “CAT” scans (computed axial tomography), to name a few.

 

·       A set of linear equations is of the form:

 

 

·       We need to determine what a, b, c & d will satisfy all four equations.

 

·       To do this in MATLAB, we need to store all the coefficients on the left hand side in a matrix (a 2D array) like so:

 

A=[2  3  1    6

      4 -2 .5    1

      6  1   0  -1

               2  2  -1   3]

 

·       Now we store the solution column in a column vector (a 1D array):

 

B=[4

      2

     -1

      4]

 

·       The solution is found using

 

x=inv(A)*B

 

·       You can also use

 

x=A^(-1)*B

 

x=A\B  % this is the most computationally efficient and stable way

 

·       The variables we want are stored in the column vector x

 

a=x(1)

b=x(2)

c=x(3)

d=x(4)

 

·       Sometimes a solution will not exist or multiple solutions will exist and the method described above will not work (inv(A) will yield an error).  In order for a single solution to exist, we need exactly N linearly independent equations for N unknowns.  An independent equation is one which cannot be formed by adding (in some combination) the other equations.

 

·       Sometimes you might not know if you have linearly independent equations just by looking at them.  The command rank( [your coefficient matrix] ) will tell how many linearly independent equations you have

 

rank(A)

 

 

Application in medical imaging

 

·       Let’s look at the tomography example.   Tomography involves the solving of a special set of linear equations.  Check out the website link below and then we will look at some of the mathematics behind Computed Tomography (CT):

 

http://www.imaginis.net/ct-scan/how_ct.asp

 

 

·       Here is a CT reconstructed image “slice” of a brain.  Note that the light areas are made up of more dense tissue than the dark regions.  Thus, the bright regions will dampen x-rays more so than the dark regions. 

 

 

 

·       Each small square above (picture element or pixel) corresponds to an initially unknown variable (a density value) that we need to solve for.  Our equations come from measuring the strength of x-rays passing through a brain at different angles.

 

·       For simplicity, imagine a very simple brain slice with only 4 distinct density regions (you may know someone like this). 

 

·       We send a thin sheet of x-rays through at different angles and measure the x-ray intensity getting through.  Based on the physics, it turns out that the measured intensity is the negative exponential of the sum of densities along the path that the sheet of x-rays takes.

 

·       Note that when the sum of density is zero (air), the measured intensity is 1 (full strength).  When the sum densities are infinite (maximum density) the measured intensity is 0.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


·       The 4 measurements (m1-m4) lead to 4 linear equations (take the natural log of both sides):

 

 

·       These are rearranged in a standard form

 

 

·       Let’s see if we have enough linearly independent equations…

 

A=[1 1 0 0

       0 0 1 1

       0 1 0 1

       1 0 1 0];

 

rank(A)

 

·       Problem!  Equation 4 can be made by adding Equations 1 and 2, and then subtracting Equation 3.  Note: adding equations means to add each corresponding term on the left and right hand sides. 

 

·       Because our rank is 3 and we have 4 unknowns, we need one measurement to be different somehow.  How about this…

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


·       We now get these linear equations:

 

 

·       Let’s see if we have enough linearly independent equations now?

 

A=[1  1  0 0

       0  0  1 1

       1 .5 .5 1

       1  0 1  0];

 

rank(A)

 

·       Let’s say we measure the following relative x-ray intensities

 

m1 = .3679

m2 = .0907

m3 = .0954

m4 = .1827

 

B=[-log(m1)

       -log(m2)

       -log(m3)

       -log(m4)];

 

·       Now we can calculate the relative densities of the four brain quadrants

 

x=inv(A)*B

 

a=x(1)

b=x(2)

c=x(3)

d=x(4)

 

·       Once the densities are known, we can form an image

 

scan=[a b

            c d]

 

imagesc(scan);     % display a scaled version of our densities

colormap(gray(256))

title('Our Computed Tomography Image');

 

·       You can now see the internal density structure in a way that a regular x-ray cannot!  Obviously, to be useful, we need to solve for more variables (more pixels) and thus we will need more detectors and more angles to give more equations.  That is how a real system works.

 

·       Just as a reality check on our answer, consider this: the blocks a & b are the least dense (darkest) and that gives rise to a largest measurement m1, because the x-rays are not passing through dense tissue.  The blocks c & d are the densest (brightest) and this means that m2 will be small because the x-rays passed through both of these dense blocks. 

 

·       In a mathematics course called “Linear Algebra,” the subject of solving multiple linear equations is treated in much greater detail.  You will also find more details in Chapter 5 (pages 193-235). 

 

·       Issues like what to do when you have too few equations, and can’t get more the way we did here by changing the experiment, are covered (the underdetermined case p. 213).  Another important case to be aware of is the case when you have too many equations.   Extra equations are helpful when noise is present in your measurements (from which your equations are derived).  In that case, none of your equations are perfectly correct, but all provide some useful information (overdetermined case p. 227).