d i
ei ady
b

omputing

Antti Rasila, Susanna Liesipohja, Juha Kuortti

May 16, 2011

1 Introdution to MATLAB 2

1.1 Basis . . . 2

1.1.1 Output. . . 4

1.2 Vetors and matries . . . 6

1.2.1 Random numbers . . . 9

1.3 Branh and loopstrutures . . . 9

1.4 Dening funtions. . . 13

1.5 Polynomials . . . 14

1.6 Plottingand drawing . . . 15

1.6.1 Plotting 3D graphis . . . 16

1.7 Useful links . . . 16

2 Linear algebra 17 2.1 Linear equations . . . 17

2.2 Matries and vetors inMATLAB . . . 19

2.2.1 Solving linear equationsin MATLAB . . . 27

2.3 Gaussianelimination . . . 28

2.4 Matrixdeompositions . . . 30

2.4.1 LU-fatorization. . . 30

2.4.2 Cholesky-deomposition . . . 34

2.4.3 QR-deomposition . . . 36

2.4.4 Singularvalue deomposition . . . 40

2.5 Linear least squares . . . 43

2.5.1 Least squares and MATLAB . . . 48

2.6 Symboli linearalgebra inMATLAB . . . 50

3.1 Polynomial interpolation . . . 53

3.1.1 Lagrangeinterpolation . . . 53

3.1.2 Determiningoeients . . . 55

3.2 Runge's phenomenon . . . 56

3.3 Pieewise linear interpolation . . . 57

3.4 Splines . . . 59

3.4.1 Cubi spline . . . 60

3.5 Additionalmethods forinterpolation inMATLAB . . . 63

4 Numerial dierential and integral alulus 67 4.1 Numerialderivation . . . 67

4.1.1 Estimatingderivative with polynomial . . . 72

4.2 Jaobianmatrix . . . 73

4.3 Numerialderivation onomplex plane . . . 75

4.4 Numerialintegration. . . 77

4.4.1 Trapezoid rule . . . 77

4.4.2 Simpson's rule . . . 79

4.4.3 Numerialintegration inMATLAB . . . 81

4.4.4 Numerialintegration onomplex plane . . . 83

4.4.5 More advaned integration methods . . . 86

4.5 Symboli dierentialand integralalulus . . . 93

5 Nonlinear equations 98 5.1 Root ndingalgorithms . . . 99

5.1.1 Braketing . . . 99

5.1.2 Fixed pointiteration . . . 100

5.1.3 Seant method . . . 103

5.1.4 False positionmethod . . . 105

5.1.5 Newton's method . . . 106

5.1.6 Brent's method . . . 110

5.1.7 Roots of polynomials . . . 114

5.1.8 Root ndingin MATLAB . . . 117

5.2 Minimizationalgorithms . . . 118

5.2.1 Golden setionsearh . . . 119

5.2.2 Brent's method . . . 121

5.2.3 Searh methodsfor multivariablefuntions . . . 124

5.2.4 Searhing minimum in MATLAB . . . 128

6.1 Numerialsolutions toODE's . . . 132

6.1.1 Euler's method . . . 133

6.1.2 Runge-Kutta methods . . . 136

6.2 Solving ODE's inMATLAB . . . 139

6.3 Boundary value problems . . . 141

6.4 Partial dierentialequations . . . 142

6.4.1 Waveequation . . . 143

6.4.2 Heat equation . . . 146

Introdution to MATLAB

MATLABisaninterativeomputingenvironmentfordoingnumerialom-

putations with vetors and matries. It is suitable for a variety of dierent

tasks involving sienti and tehnial omputations. MATLAB was re-

ated by Cleve Moler in the 1970's as a small program for teahing matrix

alulations. It was reated by using Fortran. LINPACK and EISPACK

library routines were used internally for omputations. MATLAB was ini-

tially a shareware program, and it quikly spread to other universities. In

1984, Cleve Moler, Jak Little and Steve Bangert founded MathWorks and

ommersialized MATLAB. The urrent versions are ommerial produts,

written inC. Several extensions have been addedto the originalMATLAB.

1.1 Basis

As mentioned above,MATLAB an be used interatively asa sort of alu-

lator. Inthis ase,the ommandswillbewrittendiretlyintothe MATLAB

prompt >>. One an exit the promptby writingquit orexit.

The program an also be written into a le ending with .m, for example

myfile.m, and run from the prompt by typing the lename, in this ase

myfile.

The MATLAB ommands are organized into dierent topis. Typing help

willgivealistof allthetopis andtyping help [topi℄willgivealistof all

the ommandsgrouped under that topi. Typing help [ommand℄ willgive

a short desription of the spei ommand.

TheMATLABommandsissuedandthe resultsobtainedan besaved using

>> diary test.dry

>> a=1; b=0;

>> a+b

ans =

1

>> diary off

>> type test.dry

a=1; b=0;

a+b

ans =

1

diary off

Calulations are done inoating-point preision(approximately16 digitsin

thedeimalsystem). Theoutputanbehangedusingtheformatommand,

but it will not hange the preision of the alulations. The default output

preision isshort, whihis of 5-digit preision. Forexample

>> pi

ans =

3.1416

>> format long

>> pi

ans =

3.141592653589793

A variable is given a value with the = operator. The most ommonly used

variableisans, italwaysontainsthe resultofthe previousommand. More

preisely, if a ommand does not assign a value to a named variable then it

is stored tothe variableans.

Some names are reserved for ertain onstants,suh aspi (for

π

^{),and both}

i and j represent the imaginary unit. Other reserved names are, among

others,realmax,realmin,eps,InfandNaN.Dierentvaluesanbeassigned

to these onstants, but they will revert bak to the default values after re-

starting the programor using the ommand lear.

A short example:

>> 1/0

ans =

ans =

NaN

>> NaN=5

NaN =

5

>> lear

>> NaN

ans =

NaN

1.1.1 Output

A simple way to output data is to display a variable. This an be aom-

plished by giving its name (without a semiolon) in interative mode. Al-

ternatively you an use the disp funtion, whih shows values without the

variablename, asin:

x=42;

>> disp(x)

42

For a fanier output, MATLAB has various funtions for reating strings

from numbers, formattingdata, et... One suh is fprintf, whih an also

beused for printing intoa le. The syntaxfor this is:

fprintf([fileId℄,[format℄,[ inpu t values℄)

If [fileId℄ is omitted, the funtion will print diretly onto the sreen.

[fileId℄referstotheleidentierreturnedwhenopeningtheleforwriting

with fopen. For example, fileId=fopen('myfile.txt' ,'w' ) would open

myfile.txtforwriting. Theommandflose(fileId)wouldlosethele.

[format℄ is a string in single quotation marks that desribes the format of

the outputelds. It an inludeombinationsof the following:

•

^{Aperentsignfollowedbyaonversionharater,suhas%sforstrings}

and %d for an integer. Floating-pointnumbers an beprinted with %f

for xed notation and %e for exponentialnotation.

•

^{Fieldwidthandpreision. Forexample,%6.2fwouldrefertoaoating-}

pointnumber of eld width 6and preision 2.

•

^{Flags, suh as- for}left-justied and+ for printing asign harater (

+

or

−

^{). Forexample, %+-d would print a signed integer justied to the}

left.

•

^{Literal textto print.}

•

^{Esape haraters, suh as \n for a new line\t for tab and %% for the}

perent sign.

Beloware someexamples on the use of fprintf.

a=5; s='Hello world';

>> fprintf('%d is an integer and %s is a string\n',a,s)';

5 is an integer and Hello world is a string

>> fprintf('Now %+d is a signed integer\n',a)

Now +5 is a signed integer

b=1.23456789; =0.0015;

>> fprintf('Printing with preision 2: %.2f\n',b)

Printing with preision 2: 1.23

>> fprintf('\t or with width 20: %20f\n',b)

or with width 20: 1.234568

>> fprintf('Printing as %f and as %e\n',,)

Printing as 0.001500 and as 1.500000e-03

For printing intoa le, one an do the following:

>> fid=fopen('output.txt','w');

>> fprintf(fid, '%s\n',s);

>> flose(fid);

Nowthe sentene Hello world (and a row-hange)an be found inthe le

output.txt.

In MATLAB, the basi data struture is matrix. The most eient way of

programming MATLAB is to treat every variable as a vetor or a matrix.

Assigning vetor values an be done inthe followong ways:

>>x = 1:1:4; % expression a:h:b produes a vetor with

% numbers from a to b with interval h. If

% no h is provided, 1 is assumed, eg. 1:10

>>y = [0 1 0 1℄; % Vetor values an be

% given individually also.

Vetor dimensions have to be taken into aount when performing arith-

metis. The produt x*y is not dened for two

n

^{-vetors, but the pairwise}

operations x.*yand x+yare:

>>x.*y

ans =

0 2 0 4

>>x.+y

ans =

1 3 3 5

In the ase of vetors, the produt is dened as if they are

n × 1

^-matries:

heneweneedtotransformonevetorfromarowvetortoaolumnvetor.

We dothis with the transpose operator '.

>>x

x =

1 2 3 4

>>x'

x =

1

2

3

4

>>x'*y

ans =

0 1 0 1

0 2 0 2

0 3 0 3

0 4 0 4

ans =

6

If your vetor (or matrix) ontains omplex numbers, you need to take into

aount that the transpose operator will also hange a omplex number to

its omplement, i.e., if

z = a + bi

^then

z ¯ = a − bi

^.

>> xi=[2+i 2 -i 4℄;

>> xi'

ans =

2.0000 - 1.0000i

2.0000

0 + 1.0000i

4.0000

The power operatoris ^, and again, itonly works elementwise:

>>x.^y

ans =

1 2 1 4

Elementaryfuntions are alsoavailable forvetors:

>> sin(x)

ans =

0.8415 0.9093 0.1411 -0.7568

>>exp(y)

ans =

1.0000 2.7183 1.0000 2.7183

You an dene a matrix just as you dened a vetor: to indiate a row

hange,use ;

>>A = [1 2 ; 3 4℄

A =

1 2

3 4

>>b= [5; 6℄; % b must be a row vetor

% You an now obtain inverse of A and multiply

% b with it

-2.0000 1.0000

1.5000 -0.5000

>>x = iA*b

x =

-4.0000

4.5000

% It is generally faster and easier to use MATLABs

% built-in linear solver operator \

>> x = A\b

x =

-4.0000

4.5000

Some useful matrix ommands are also: eye (produes an identity matrix),

zeros (produes a matrix of all zeros) and ones (produes a matrix of all

ones). It is alsopossible to selet spei elements, rows oh olumns from

a matrix. The ommand for this is A[i,j℄, where A represents a matrix, i

the row of that matrix and j the olumn. Here, i and j an be salars or

vetors.

% We reate a 3x3-matrix of all ones

>> A=ones(3)

A =

1 1 1

1 1 1

1 1 1

% To pik a speifi element from matrix A, use A(i,j).

% To pik a whole row (or olumn), replae j (or i)

% with :

>> A(1,:)

ans =

1 1 1

>> A(2,:)=[2, 3, 4℄

A =

1 1 1

2 3 4

1 1 1

A =

1 1 1

2 3 4

1 42 1

1.2.1 Random numbers

Randomnumbers anbegeneratedby usingtheommandsrandandrandn.

The ommand rand(m,n) will produe an

m × n

^{-matrix of uniformly dis-}

tributed randomnumbers on

(0, 1)

^{and randn(m,n)willprodueamatrixof}

normally distributed random numbers with mean

0

^{and standard deviation}

1

^.

dist = zeros(6,1);

for j=0:99

k = round(5*rand(1)+1);

dist(k) = dist(k)+1;

end

disp(dist)

Output:

1: 14

2: 15

3: 13

4: 18

5: 18

6: 22

1.3 Branh and loop strutures

The branh and loopstrutures availablein matlabare: for, while, if and

swith. Themainprinipleisthatyoushouldonlyusetheseasalastresort.

If possible, you should use eient vetor operations instead.

The syntaxof the forstatement is:

for [variable℄=[vetor℄

...

Ifonewantstorepeattheloop

k

^{times,itishandytousethevetorstatement}

1:k, whih produes alist of numbers

1, 2, . . . , k

^.

# Example: 2nd powers of positive integers

for x = 1:4

xx = x*x

fprintf('%d * %d = %d',x,xx)

end

Output:

1*1 = 1

2*2 = 4

3*3 = 9

4*4 = 16

The syntaxof the if statement inMATLAB is:

if [ondition℄

...

elseif [ondition℄

...

else

...

end

The elseif and else branhes may be omitted. The ommands in the if

branh are exeuted if the ondition is satised, if not then the onditions

in the elseif branhes are evaluated. If none of the onditions given is

satised, the ommands inthe else branhare exeuted.

The most ommononditions used are the omparisonoperations <, <=, ==,

=, >= and >. Note that for equality, the expression == is used in order to

avoidonfusion withthe value assignmentoperator= 1

. The expression = is

used for inequality.

A while statement is used when one wishes to repeat the loop until some

ondition is no longer satised. This struture is very useful when reading

input from ale orfromthe user.

The syntaxof the while statement is:

while [ondition℄

...

end

1

Thisissigniant,asine.g. theClanguage,theonditionif(x=1)...isalwaystrue.

to invalidate the ondition when the desired number of loops is reahed.

x = 5;

guess = 0;

while guess ~= x

guess = input('Guess a number:');

if (abs(guess - x)>10)

disp('Your guess is very wrong')

end

end

Output:

Guess a number:6

Guess a number:100

Your guess is very wrong

Guess a number:5

The syntaxof the swith statementis:

swith [swith expression℄

ase [ase expression 1℄

...

ase [ase expression 2℄

...

...

otherwise

...

end

The statements assoiated with a ertain ase will be exeuted when the

swith expression equals the ase expression in question.

% Color evaluation

olor='aqua';

swith olor

ase {'red','pink','rose'} % multiple ase expressions

disp('The olor is red.')

ase {'blue','turquoise','aqua'}

disp('The olor is blue.')

ase 'yellow'

disp('The olor is yellow.')

end

Remark. One should avoid omparing non-integers with the == operator.

For example, pi==3.14159265... is atually false. The MATLAB pi is

onlyalulatedtoaspeilength,andthus, doesnotatuallyequal

π

^{. The}

followingprogram willdemonstratethis fat:

% Desired auray of approximation

tol=10^-4;

mypi=1;

while mypi~=0

mypi=input('Guess the value of pi (0 exits):');

if mypi==pi

disp('Comparison to MATLAB pi is true')

else

disp('Comparison to MATLAB pi is false')

end

if (abs(pi-mypi)<tol)

disp('Close to pi!')

elseif (abs(pi-mypi)>1)

disp('Far from pi!')

else

disp('Not lose enough to pi!')

end

end

d Output:

Guess the value of pi (0 exits):3.141592653589793238

Comparison to MATLAB pi is true

Close to pi!

Guess the value of pi (0 exits):3.141592653589793258

Comparison to MATLAB pi is true

Close to pi!

Guess the value of pi (0 exits):3.14

Comparison to MATLAB pi is false

Not lose enough to pi!

The rst guess isan aurate approximation of

π

^{, but the seond one is not}

(the seond-to-last digit is wrong). However, the omparison to MATLAB

pi isorret inboth ases.

A funtion an be dened with the funtion statement. The syntax of this

statement is:

funtion [output℄=[funtion name℄([input℄)

...

This funtion should be saved inan m-le and the nameof the lemust be

the same as the funtion name. For example, the funtion below should be

saved assolve2.m.

funtion x = solve2(a,b,)

D= b^2 - 4*a*;

% Floating point number should not be diretly

% ompared to zero

if(abs(a)<1e-6)

disp('Error')

return

else if(abs(D)<1e-6)

x = -b/2*a;

return

else

x(1) = -b + sqrt(D)/2*a;

x(2) = -b - sqrt(D)/2*a;

end

end

Output:

>> solve2(1,0,0)

ans =

0

>> solve2(1,0,1)

ans =

0 + 1.0000i 0 - 1.0000i

>> solve2(1,0,-1)

ans =

1 -1

Theabovefuntionsolve2solvesthe rootsofagiven seondorderequation.

The input parameters given for the funtion are three numbers

a

^,

b

^and

c

^,

example,theequation

x ² = 0

^{issolved(onerootat}

0

^{);intheseondase,the}

equationis

x ² + 1 = 0

^{(onlyimaginaryroots)andinthethirdase}

x ² − 1 = 0

(tworoots

± 1

^).

Forsimpler funtions,itmay beeasier todenethe funtions diretly into

the program. Thisanbedonewiththe inlineommandor,morereently,

the funtionhandle .

>> f=inline('exp(x.^2)','x')

f =

Inline funtion:

f(x) = exp(x.^2)

>> g=(x) x.^2

g =

(x)x.^2

1.5 Polynomials

In MATLAB, a polynomial is represented by a vetor whih onsists of its

oeients. To reate a polynomial one an simplyenter eah oeient of

the polynomial into the vetor in desending order. For instane, onsider

the following polynomial:

p(x) = 2x ⁴ − x ² + 5x + 17

To givethis in MATLAB, justwrite the vetor

>> p=[2 0 -1 5 17℄;

One may nd the value of a polynomial using the polyval funtion. For

example, tond the value of the above polynomial at

x = 2

^,

>> polyval(p,2)

ans =

55

Therootsofapolynomialanbeobtainedwithroots([your polynomial℄).

For example,the rootsof the polynomialabove are

>> r=roots(p)

r =

1.2663 + 1.3591i

1.2663 - 1.3591i

-1.2663 - 0.9273i

Ifoneknows theroots already,theoeientsan befound usingtheinverse

funtionpoly. Twopolynomialsanbemultipliedbyusingonv([poly1℄,[poly2℄)

and dividing an bedone in asimilar way with the deonv funtion.

1.6 Plotting and drawing

Curves an be drawn with the ommand plot. For example, to plot a sine

urve, one an do the following:

>> x=0:.1:2*pi;

>> plot(x,sin(x))

It is possible to plot several urves at one. The appearane of the urves

an behanged. For example,the ommand

>> plot(x,sin(x),'r',x,os(x),'.b')

plots the sine urvein red and the osineurve as blue dots.

0 1 2 3 4 5 6 7

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

Onemaylabeltheaxeswiththeommandsxlabel('[labelname℄'),forthe

legend([urve1℄,[urve2℄,[ urve 3℄,[ ... et.℄).

1.6.1 Plotting 3D graphis

Spatial urves given in parameters an easily be plotted with the funtion

plot3 simplyby adding

z

-oordinates.

Surfae plotting an be done with the funtion surf. But rst, one should

generate the appropriateX and Y arrays using the funtionmeshgrid.

Inthe examplebelow, whereweare plottingthe funtion

f (x, y) = xe ^{− x 2 − y 2}

^,

X and Y represent the plane and Z represents the height.

>> x=-2:.1:2; y=x;

>> [X,Y℄=meshgrid(x,y);

>> Z=X.*exp(-X.^2-Y.^2);

>> surf(X,Y,Z)

Output:

−2

−1 0

1

2

−2

−1 0 1 2

−0.5 0 0.5

1.7 Useful links

Linear algebra

2.1 Linear equations

An equation with variables

x 1 . . . x n

^{that an be writtenin the form}

a 1 x 1 + a 2 x 2 + · · · + a n x n = b, n ≥ 1

isalledlinearequation. Theoeients

a 1 . . . a n

^and

b

^{anberealoromplex}

numbers.

A system of linear equations is a olletion of one or more linear equations

involvingthe same variables. Using matrix algebra,the linear system



 

 



 

 



a ₁₁ x ₁ + a ₁₂ x ₂ + . . . + a _1n x _n = b ₁ , a ₂₁ x ₁ + a ₂₂ x ₂ + . . . + a _2n x _n = b ₂ ,

.

.

.

.

.

.

.

.

.

a m1 x 1 + a m2 x 2 + . . . + a mn x n = b m ,

an bewritten inthe form

Ax = b; A =











a 11 . . . a 1n

a 21 . . . a 2n

.

.

.

.

.

.

a _m1 . . . a _mn









 , x =









 x 1

x 2

.

.

.

v _n









 , b =









 b 1

b 2

.

.

.

b _m









 .

The matrix

A

^{is a}

m × n

^{matrix, the vetor}

x

^{avetor with}

n

^omponents

Let

A ∈ C ^{n × n}

. A is said to be invertible, if there exists suh

B ∈ C ^{n × n}

,

alledtheinverse ofmatrix

A

^{,andisdenoted}

A ^{− 1}

^{. If}

A ^{− 1}

^exists,

A

^isalled

invertible. Usingthis denition, we get followingtheorem.

Theorem 2.1. The linear system of equation

Ax = b, A ∈ C ⁿ × n , x ∈ C ⁿ , b ∈ C ⁿ

has a single solution only if

A

^is invertible. The solution is

x = A ^{− 1} b

^.

Generally, system of linear equations an have an exat solution only if it

has exatlyas many linearilyindependent equationsas it has unknowns. In

this situationthe system of linear equations



 

 



 

 



a 11 x 1 + a 12 x 2 + . . . + a 1n x n = b 1 , a 21 x 1 + a 22 x 2 + . . . + a 2n x n = b 2 ,

.

.

.

.

.

.

.

.

.

a _n1 x ₁ + a _n2 x ₂ + . . . + a _nn x _n = b _n ,

translates into a

n × n

^{matrix and}

n

^{vetors. However, numberof equations}

and unknowns do not always oinide. In this ase we get asystem



 

 



 

 



a 11 x 1 + a 12 x 2 + . . . + a 1n x n = b 1 , a 21 x 1 + a 22 x 2 + . . . + a 2n x n = b 2 ,

.

.

.

.

.

.

.

.

.

a m1 x 1 + a m2 x 2 + . . . + a mn x n = b m , .

This system an be writtenas amatrix equation

Ax = b ,

where

A

^{is a}

m × n

^matrix,

x

^an

n

^{-vetor, and}

b

^a

m

^-vetor.

If

m < n

^{, that is, if there are fewer equations than there are unknowns,}

system is alled underdetermined. Solving an underdetermined system of

equations will not usually produe an exat solution, but the solution will

have degrees of freedomdepending onthe oeient matrix: the number of

whih is determined by how many unknowns remain in the solution vetor

x

^{. The solution an be} interpreted as a spae where the objets dened by the equationsinterset.

Ax = b =





1 3 3 2

2 6 9 5

− 1 − 3 3 0











 x 1

x 2

x 3

x 4









=



 1 5 5



 .

Using elementaryrow operationsweget the solution:

x =









− 2 − 3x ₂ − x ₄ x 2

1 − ^x 3 ⁴

x 4







 .

The vetor

x

^{is a solutionfor the equation}

Ax = b

^{with arbitrary values of}

x 1

^and

x 4

^{, thus givinganinnitely many solutions.}

If in system ??

m > n

^{, the system is alled} overdetermined; that is, there are more equations in the system than there are unknowns. Computing an

exat solution to a overdetermined system may be possible. However, the

more there are onstraints (equations), the less likelyit isthat they allhold

for a spei point. Thus solving overdetermined systems usually inludes

searhing a best possible solution: asolution that does not neessarily hold

for all the equations in the system, or any of them, but it almost holds for

all of them. This is usually ahieved by linear least squares, whih will be

disussed indepth inlater hapters.

Later we willintrodue some methods for nding the inverse of a matrix in

MATLAB.

2.2 Matries and vetors in MATLAB

ThebasidatatypeinMATLABisarealvaluedmatrix,anddefaultassump-

tion for every operation isthat the operands are matries. Some operations

are, howeverdened alsoelementwise, soastomakeertainoperationseasy

and eient, and it requires a ertain amount of alertness to avoid any ob-

vious pitfalls.

Thebasimultipliation,denotedby*, ismatrixmultipliation. Itisdened

for matries

A

^and

B

^{, where}

A ∈ C ^{m × n}

and

B ∈ C ^{n × k}

. It also works on a

obtain an elementwise produt with operator.*. The elementwise operator

willprodue the Hadamard produt of two matriesof same dimensions.

% Define two square matries and two

% non square matries

>> A = [3 2 3; 3 4 3; 4 5 1℄

A =

3 2 3

3 4 3

4 5 1

>> B = [1 2 4; 1 4 6; 1 7 7℄;

>> D = [1 2; 4 3 ; 7 6℄

D =

1 2

4 3

7 6

>> K=[3 4 5; 5 6 7℄

K =

3 4 5

5 6 7

% Multipliation of two square

% matries works

>> A*B

ans =

8 35 45

10 43 57

10 35 53

% Hadamard produt of two square

% matries

>> A.*B

ans =

3 4 12

3 16 18

4 35 7

% Produt of 3x2 and 2x3 matries

>> D*K

ans =

13 16 19

27 34 41

% Elementwise produt doesn't work

>> D.*K

??? Error using ==> times

Matrix dimensions must agree.

% Multiply B by four

>> 4*B

ans =

4 8 16

4 16 24

4 28 28

Same applies tothe power operator: the operator literallymeans that the

rst operand is multiplied by itself as many times as the seond operand

orders. Ifthe rstoperandnotasquarematrix, theoperationisnot dened.

Thusthethepoweroperatorshouldbeusedonlyasanelementwiseoperation:

.

% Examples of the power operator:

% First on real number

>> 3^5

ans =

243

% then on a square matrix

% Note that this is a defined

% operation beause A*A is a

% defined operation.

>> A^2

ans =

27 29 18

33 37 24

31 33 28

% We try then the elementwise

% power operator. Notie the

% differene with the regular

% power operator.

>> A.^2

ans =

9 4 9

% Power operator doesn't work on a non

% square matrix beause D*D is not defined

>> D^4

??? Error using ==> mpower

Matrix must be square.

% However, an elementwise operator is defined:

>> D.^4

ans =

1 16

256 81

2401 1296

The usual division sign - /, should be used only on matries with a single

value. In ase of single value, it works as one would expet: it performs

a division. However, if given matrix values, the values it produes are not

what one would expet, and obtainable inmuh more intuitiveway through

the bakslash-operator, whih we will disuss later. Those interested in us-

ing it should familiarize themselves with the mldivide manual page. The

elementwise version of division-operator is ./, whih is useful on a number

of oasions. Beause sometimes both elementwise operation, and matrix

operation an be invoked, autionis required.

Additionandsubtrationare elementwiseoperations:

A + B

^{isthe standard}

matrix addition, whih requires that both

A

^and

B

^{have the same dimen-}

sions. The addition and subtration operators have alsobeen overloaded to

inlude operations like

2 + A

^{. This operation is dened as add 2 to every}

elementof

A

^{. There,however, is not an}elementwise operation,that would allowone to addtwomatrieshaving the same numberof elements,but dif-

ferent dimensions, together. All of the above holds for the subtration as

well.

% Examples of addition

% Sum of two matries of equal sizes is ok

>> A+B

ans =

4 4 7

4 8 9

5 12 8

% So is adding 2 to every element of A

ans =

5 4 5

5 6 5

6 7 3

% This doesn't work beause K and L have

% different dimensions

>> K = [1 2 3℄

K =

1 2 3

>> L = [1;2;3℄

L =

1

2

3

>> K+L

??? Error using ==> plus

Matrix dimensions must agree.

Most of MATLAB's built-in funtions, like exp, sin and os are dened

elementwise.

% Define an even spaed real valued vetor H

>> H = 1:0.5:3

H =

1.0000 1.5000 2.0000 2.5000 3.0000

% Take an sin of eah element of the vetor

>> sin(H)

ans =

0.8415 0.9975 0.9093 0.5985 0.1411

Another topi that will require some attention is the matrix and vetor di-

mensions. As mentioned, almost all the operations are dependent on the

dimensions of the operands. Oftentimes, like when rafting a funtion, one

does not wish tox the matrix dimension,but dynamially adaptto the di-

mensions. Thewaytodothisistousefuntionslength andsize. Funtion

length is primarily meant for work with vetors, and it returns the largest

dimension of argument. For example length(ones(4,2)) would return 4.

The funtion size returns a vetor ontaining all dimensions. It is more

versatile than length, but to work,it requires anassignment. Forexample,

a

>> rows = dims(1);

Another operation that is frequently needed in order to handle the dimen-

sionsisthetranspose. MATLABdefaultsthetransposetoonjugateversion,

working astranspose onrealmatries,butreturningthe onjugatetranspose

on omplex matries. The onjugate transpose operator is the '. If one

wishes to obtain a non-onjugate transpose, a funtion transpose is avail-

able. For work on more omplex strutures than two-dimensional arrays,

MATLAB providesthe funtionpermute.

% Create a omplex matrix C

% Reall that i is overloaded

% to at as a omplex oeffiient

>> C = A+B*i

C =

3.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 4.0000i

3.0000 + 1.0000i 4.0000 + 4.0000i 3.0000 + 6.0000i

4.0000 + 1.0000i 5.0000 + 7.0000i 1.0000 + 7.0000i

% Notie the onjugate or hermitian transpose,

>> C'

ans =

3.0000 - 1.0000i 3.0000 - 1.0000i 4.0000 - 1.0000i

2.0000 - 2.0000i 4.0000 - 4.0000i 5.0000 - 7.0000i

3.0000 - 4.0000i 3.0000 - 6.0000i 1.0000 - 7.0000i

% Should you ever need it, a non hermitian transpose

% is also available.

>> transpose(C)

ans =

3.0000 + 1.0000i 3.0000 + 1.0000i 4.0000 + 1.0000i

2.0000 + 2.0000i 4.0000 + 4.0000i 5.0000 + 7.0000i

3.0000 + 4.0000i 3.0000 + 6.0000i 1.0000 + 7.0000i

% Transpose of a real valued vetor

>> K'

ans =

1

2

3

The matriesan be indexedwith two numbers, asusual, the rst being the

row-index, the seond being the olumnindex, and indexes starting from1.

way to index matries. Matries an be indexed with a single number, the

index runningdown olumn wise. That is, A(3) = A(3,1). While one may

do this, for the sake of larity, it is highly disouraged. The reason this

option is available is due to the properties of omputer arhiteture and C,

the language that MATLAB iswritten with.

MATLAB allows aessing entire rows, olumns, and submatries of any

matrix. Thisisahievedwiththe rangeoperator:. Ifnotgivenanyrange,it

defaults towholeroworolumn,for example: the ommand A(:,1)returns

therstolumnof

A

^{,whileA(2,:) wouldreturntheentireseondrowof}

A

^.

Instead of seleting the entire row or olumn, one an selet only a part of

it by givingthe rangeoperator parameters: A(1:5,1) wouldreturn the rst

ve elements the rst olumn of

A

^{. Seletion of submatries follows suit:}

insteadof givingonerange, wegivetwo: A(2:3,3:4) would returnamatrix

that wouldontain

A

^'selements

a 2,3 , a 2,4 , a 3,3

^and

a 3,4

^{. Thisseletionanbe}

extended further: seletion index an be any olletion of positive integers,

and the seletionstill works, as long as they are within index bounds of A,

for example seletionA([1 3 5℄,2) returnselements

a 12 , a 32

^and

a 52

^.

Seletion methods are not limited to numerial indexing; it is also possible

to invoke so alled logial indexing. Logial indexing is ahieved through

reatingalogial array,and givingitasaindex. Logial arrays are returned

by logial operators, & ,| and , relational operators, suh as ==, =,>

and <, as well as any logial funtions, suh as any, isinf and isequal.

Using these operators and logial indexing, we an, for example, selet all

the positive elementsof a matrix.

% Define a large enoug a matrix

>> A = [-2 3 2 4 -4 0; -3 -4 -5 -11 2 4

3 -5 3 2 3 4; 1 -3 2 -4 5 -6; 1 2 3 -4 6 5℄

A =

-2 3 2 4 -4 0

-3 -4 -5 -11 2 4

3 -5 3 2 3 4

1 -3 2 -4 5 -6

1 2 3 -4 6 5

% Selet the third row of the matrix

3 -5 3 2 3 4

% The alternate indexin way: A(12) is the

% the same as A(rem(12,5),mod(12,5)+1).

% While there are situations it an be

% more effiient than the usual way, readability

% suffers.

>> A(12)

ans =

-5

% Seleting submatries is quite similar to

% single elements or rows and olumns: just

% give to ranges

% Here we have seleted rows 2 3 4 and 5, and

% olumns 3 4 5 and 6.

>> A(2:5,3:6)

ans =

-5 -11 2 4

3 2 3 4

2 -4 5 -6

3 -4 6 5

% Finally a look into the logial seletion routines:

% selet all the elemnts of A less than -4

>> A(A<-4)

ans =

-5

-5

-11

-6

% A more ompliated logial ondition: selet

% elements of A smaller than 0 but greater than

% -5.

% Note that this requires the use of a elementwise

% logial operator &, whih is defined for the use

% with logial matries and vetors.

>> A((A<0)&(A>-5))

ans =

-2

-3

-4

-4

-4

-4

There are several matries, that ome up often in linear algebra, most no-

tablebeingtheunitmatrix. MostoftheseareprovidedinMATLAB'smatrix

library,whihgenerates themaording togiven parameters. The ommand

eye produesthe unit matrix of given dimensions. The ommand ones pro-

dues a matrixomposed entirelyof ones, and the ommand zeros, aord-

ingly,produes amatrix made up of zeros. Some of the moreexoti built-in

matries are, for example, the Hilbert matrix and the magi square. The

Hilbert matrix is produed by ommand hilb. The Hilbert matrix is de-

signedtohaveertainverypoornumerialproperties. Themagisquareisa

square matrixwithequalolumn,rowand diagonalsums,and itisprodued

by the ommand magi.

2.2.1 Solving linear equations in MATLAB

The primarytoolfor solving linearequationsinMATLAB isthe \-operator.

Tosolvealinearequationoftheform

Ax = b

^{weusetheommandx = A \b.}

The bakslash operator isvery versatile: if the matrix

A

^is overdetermined, i.e, there are more rows than there are olumns, a solution in least-square

senseisprovided. Ifthesystemisunderdetermined,itndsthebasisolution

with at most

m

^{nonzeros. Here are afew examples:}

% Example onerning the Hilbert matrix

>> A= hilb(10);

>> x = ones(10,1);

>> b = A*x;

>> sol = A\b;

>> norm(x-sol)

ans =

8.7188e-04

% The previous lsq-example with the bakslash-operator

>> A = [1 1; 2 1; 3 1 ; 5 1; 7 1; 9 1 ; 10 1℄;

>> b = [444 458 478 506 523 543 571℄;

x =

13.0798

434.1498

2.3 Gaussian elimination

Probably the most famous method for solving an

n × n

^{system of linear}

equationsis theGaussian elimination algorithm,named afterCarlFriedrih

Gauss. The idea of the algorithm is to, for eah olumn of the oeient

matrix, eliminatethe elementsbelowthe diagonalusingrowoperations,and

when anuppertriangularmatrix isahieved, we doa bakward substitution,

solving

x n

^{fromthelastequation,and}substitutingthe solutiontotheseond last equation,and thus gainingsolution to

x n − 1

^{, and soforth.}

Example2.3. SolveasystemoflinearequationsusingGaussianelimination,

when the system is:



 

 

3x 1 − x 2 + x 3 = 2,

− x 1 + 3x 2 − 2x 3 = 1, 2x 1 + 2x 2 − x 3 = − 3,

.

Eliminate all the elements below the rst element on the rst olumn: we

addthe rstrowmultipliedby

1

3

^{totheseond. Thenaddtherst multiplied}

by

− ² ₃

^{to the third row, and we get:}



 

 

3x 1 − x 2 + x 3 = 2, 0x 1 + 8x 2 − 5x 3 = 5, 0x 1 + 8x 2 + x 3 = − 13,

.

Then eliminate the seond element of the third row by adding the seond

row multipliedby

− 1

^{to it,and youget:}



 

 

3x 1 − x 2 + x 3 = 2, 0x 1 + 8x 2 − 5x 3 = 5, 0x ₁ + 0x ₂ + 6x ₃ = − 18,

.

We then obtain

x ₃ = − ¹⁸ ₆ = − 3

^{, and plae it in the equation on the seond}

row, and get

x 2 = − ⁵ ₄

^{,and nallywe get}

x 1 = ⁵ ₄

^.

The algorithmfor Gaussianeliminationin MATLAB ode is:

Listing2.1: Algorithm for Gaussianelimination

funtion x = gauss_el(A,b)

n = length(A);

% part a - elimination

for i = 1:n-1

for j = i+1:n

% alulate sale fator

m = A(j,i)/A(i,i);

% perform row operation:

% eliminate the elements below diagonal

% on olumn i

A(j,:) = A(j,:) - m*A(i,:);

b(j) = b(j) - m*b(i);

end

end

% part b - bakward substitution

x = zeros(n,1);

x(n) = b(n)/A(n,n);

for i = n-1:-1:1

x(i) = (b(i) - A(i,i+1:n)*x(i+1:n))/A(i,i);

end

Gaussianeliminationispronetonumeriinstabilitywhenworking onnearly

singular matries. The Hilbert matrix is one example of a nearly singular

matrix. Problems rise if at some part of the algorithm the absolute value

of the divisor

ˆ a kk

^{(i.e. a diagonal element after k-steps of} elimination) is verysmall. Thiseasilyleads tolossofpreisionduetothe natureof oating

point arithmeti, and auses the error to aumulate. These situations an

generallybeavoidedthroughpivotingthematrix,that is,hangingtheorder

of rows and/or olumns, and applying the same permutations both to the

solution vetor and the right-handside of the equation.

Example2.4. WenowestablishwhyGaussianeliminationwithoutpivoting

is not a stable algorithm. The funtion gauss_elim is the same as the

b

^where

A

^{is a Hilbert matrix (the Hilbert matrix is omposed as follows:}

H ij = 1/(i + j − 1)

^{) using Gaussian}elimination.

>>X = ones(13,1);

% We set up a syntheti solution to be a

% vetor omposed of ones

>>A = hilb(13);

% MATLAB provides some speial matries ready,

% Hilbert's is one of them

>>b = A*X

>>sol = gauss_elim(A,b);

% sol now holds the solution yielded by gauss_el

>> norm(sol-x)

ans =

11.0527

As is obvious, the Gaussian elimination does not provide aurate results

when dealing with matries that are badly onditioned. In numeri ases,

it is reommended to use the matrix deompositions, whih we will disuss

next.

2.4 Matrix deompositions

It is often diult to solve the equation

Ax = b

^{. Therefore in numeri}

matrix omputationwe usuallytrytopresent

A

^{asaprodutof twoormore}

matries of some simpler form. This kind of representation is alled matrix

deomposition. As we willsee, matrix deompositions will oftengive us not

only aneasier way to solvethe linearsystem, but giveus informationabout

the deomposed matrix aswell.

2.4.1 LU-fatorization

In the Gaussian elimination the matrix

A

^{is rst redued into an upper}

triangularform,from whihit iseasy toobtain solutionsthrough bak sub-

stitution. The ideabehindthe LU-fatorizationistopresent

A

^asaprodut

of two matries,

L

^and

U

^{, of whih}

U

^{is upper} triangular, and

L

^{is lower}

triangular. Wethenansolvetheequation

Ax = b

^{bysolvingtwotriangular}

( Ux = z

Lz = b ,

^{that is,}

Ax = LUx = Lz = b

TheworkingideaoftheLU-algorithmistoperformtheGaussianelimination

algorithmonmatrix

A

^{,andtakereordofthemultiplierthatwasusedtozero}

the elementsbelowthe diagonal oneaholumn. Here is aquik example:

Example 2.5.

A =





1 − 1 3 3 − 5 12 0 2 − 10



 .

Weseethat inordertoeliminatetheelements

a ₂₁

^and

a ₃₁

^{therst rowmust}

bemultipliedby3and 0respetivelybeforesubtratingfromtheseondand

third rows. Thus we get





1 − 1 3 (3) 8 3 (0) 2 − 10





where numbers inparenthesis represent the reorded multipliers. These will

form the lower triangularmatrix

L

^{. On the seond step we get}





1 − 1 3 (3) 8 3 (0) ( ¹ ₄ ) − ⁴³ ₄



 .

The diagonal elements an be inluded either in

L

^or

U

^{. The other matrix}

willhave ones onthe diagonal. Now wehave the

L

^and

U

^,

L =





1 0 0 3 1 0 0 ¹ ₄ 1



 , U =





1 − 1 3

0 8 3

0 0 − ⁴³ 4





that satisfy

LU = A .

When doingalulationswithpaperandpen, itisgenerallyeasiertouse the

so alled Doolittle algorithm. In this algorithm,the diagonalelements of

L











a 11 a 12 . . . a 1n

a ₂₁ . . . a _2n

.

.

.

a n1 . . . a nn











=











1 0 . . . 0

l ₂₁ 1

.

.

.

1

l n1 . . . l n(n − 1) 1





















u 11 u 12 . . . u 1n

0 u 22 u 23

.

.

.

0 . . . u (n − 1)n

0 0 u nn









 .

Here isan example.

Example 2.6. Let's formthe LU deomposition for the matrix

A

^{, when}

A =





6 5 12

30 18 51

− 24 − 76 − 98



 =





1 0 0

l 21 1 0 l ₃₁ l ₃₂ 1









u 11 u 12 u 13

0 u 22 u 23

0 0 u ₃₃



 .

The

3 × 3

^{matrix gives us 9 equations, eah with only one unknown. From}

the rst rowwe get

u 11 = 6, u 12 = 5, u = 13 = 12.

On the seondrow, we get



 

 

l 21 u 11 = 30 ⇔ l 21 = 5,

l 21 u 12 + 1 · u 22 = 18 ⇔ u 22 = 18 − l 21 u 12 = − 7, l 21 u 13 + 1 · u 23 = 51 ⇔ u 23 = 51 − l 21 u 13 = − 9,

and on the third



 

 

l 31 u 11 = − 24 ⇔ l 31 = − 4

l 31 u 12 + l 32 u 22 + u 23 = − 76 ⇔ l 32 = ( _u ¹ ₂₂ ) − (76 − l 31 u 12 ) = ( ₋ ¹ ₇ )( − 76 − ( − 4 · 5)) = 8 l 31 u 13 + l 32 u 23 + u 33 = − 98 ⇔ u 33 = − 98 − (l 31 u 13 ) − (l 32 u 23 ) = 22

and thus:

A = LU =





6 5 12

30 18 51

− 24 − 76 − 98



 =





1 0 0 5 1 0

− 4 8 1









6 5 12 0 − 7 − 9 0 0 22



 .

TheDoolittlealgorithmstops,ifthereappearsazeroelementonthediagonal

of

U

^{, but itis not limited toinvertible matries, in fatit an be omputed}

on matries

C ∈ C ^{m × n}

. In this ase the

L ∈ C ^{m × m}

and

U ∈ C ^{m × n}

,and the elementsbelow

u kk , k = 1 . . . m

^{, willbe zeros.}

MATLAB anompute the LU-fatorizationonanyomplex matrix

A

^with

the ommand lu. Theresult, however, is nottrue alowertriangularmatrix:

MATLAB permutes theparameter matrix

A

^{soas toahievemaximum e-}

ieny, and the

L

^{it givesis the produt of the} permutationmatrix and the atual

L

^{. To get true lower- and upper triangular matries, we get a third}

return value: the permutation matrix

P

^.

Example2.7. HereisanexampleonhowtouseLU-deompositioninMAT-

LAB.

>>A = [-1 1 4;-6 -4 0; 0 4 1℄

A =

-1 1 4

-6 -4 0

0 4 1

>>[l u ℄ = lu(A)

l =

0.1667 0.4167 1.0000

1.0000 0 0

0 1.0000 0

u =

-6.0000 -4.0000 0

0 4.0000 1.0000

0 0 3.5833

>>norm(l*u-A)

ans =

1.1102e-16

>>[l u p ℄ = lu(A)

l =

1.0000 0 0

0 1.0000 0

-6.0000 -4.0000 0

0 4.0000 1.0000

0 0 3.5833

p =

0 1 0

0 0 1

1 0 0

>>norm(l*u-p*A)

ans =

1.1102e-16

2.4.2 Cholesky-deomposition

Another matrix deomposition is the Cholesky-deomposition, named after

André-Louis Cholesky. It is not as general as LU-deomposition, but the

numberofomputationsrequiredinordertodothedeompositionissmaller.

The matries itan deompose are alsoommonin real-lifeappliations.

Denition 2.8. Matrix

A ∈ C ^{n × n}

is said to be Hermitian if it holds true that

A ^∗ = A

^{, i.e.,}

A

^{is its own onjugate transpose. This is analogialto}

the symmetry of the real matries. Note that MATLAB's '-operator gives

youthe onjugatetranspose.

Denition2.9. AHermitianmatrix

A ∈ C ⁿ × n

ispositivedenite if

h u , Au i >

0

^forall

u ∈ C ⁿ \{ 0 }

^.

A matrix

A ∈ C ^{n × n}

that ishermitianand positivedenitean bepresented with a single uppertriangularmatrix, asa produt

A = U ^∗ U.

When you have this

U

^{, you an simply use the method presented in LU-}

deomposition to solve the linear system

Ax = U ^∗ Ux = U ^∗ z = b

^{. Beause}

U ^∗

^{isa triangular matrix, this an be solved through bak} substitution.

The algorithmto produe

U

^is:

u kk = v u u t a kk −

k − 1

X

l=1

| u lk | ²

u _kj = 1 u kk

a _kj −

k − 1

X

l=1

u _lk u _lk

.

Ifthenumberunderthesquarerootisevernegative,

A

^{isnotpositivedenite}

andthealgorithmhalts,thusmakingtheCholesky-deompositionaneient

toolinstudying the positivedenity of the matrix.

Example2.10. ComputetheCholesky-fatorizationofthematrix

A

^,when,

A =





13 11 6 11 11 4 6 4 10



 .

Values of uppertriangular matrix

U

^{an be omputed withthis table:}

entry general formula value for

U

u ₁₁ √ a ₁₁ √

13

u 12 a 21 /u 11 √ 11

13

u 13 a 31 /u 11 √ 6

13

u 22

p a 22 − u ² ₂₁ q

11 − ( ¹²¹ ₁₃ ) u ₂₃ (a ₃₂ − u ₁₂ u ₁₃ )/u ₂₂ (4 − ^{√ 6} ₁₃ · ^{√ 11} ₁₃ )/ q

121 13

u 33

p a 33 − u ² ₁₃ − u ² ₂₃ q

10 − ³⁶ ₁₃ − (11 − ¹²¹ ₁₃ )

Youget anuppertriangularmatrix

U

^:

U =











√ 13 ^{√ 11} ₁₃ q

11 − ( ¹²¹ ₁₃ )

0 q

11 − ( ¹²¹ ₁₃ ) (4 − ^{√ 6} ₁₃ · ^{√ 11} ₁₃ )/ q

121 13

0 0 q

10 − ³⁶ ₁₃ − (11 − ¹²¹ ₁₃ )









 ,

having the property

A = U ^T U

^.

To obtain the Cholesky-deomposition inMATLAB, use the funtionhol:

>>A = [4 3 6; 4 7 6; 6 2 14℄

A =

4 3 6

4 7 6

6 2 14

>>A=A*A'

61 73 114

73 101 122

114 122 236

>>u = hol(A)

ans =

7.8102 9.3467 14.5962

0 3.6931 -3.9062

0 0 2.7735

>>norm(u'*u-A)

ans =

1.8201e-14

2.4.3 QR-deomposition

Any omplex square matrix

A ∈ C ^{n × n}

an be deomposed as

A = QR,

where

Q

^{is a unitary matrix, and}

R

^{is a upper triangular matrix. If}

A

^is

nonsingular,thenthefatorizationisunique,ifitisrequiredthatthediagonal

elementsof

R

^{are positive.}

Denition 2.11. A matrix

A ∈ R ^{n × n}

is orthogonal if

AA ^T = A ^T A = I

^.

A matrix

A ∈ C ^{n × n}

is unitary if

AA ^∗ = A ^∗ A = I

^.

QR-deompositionis often used tosolve problems in leas square sense. It is

the used inan algorithmfor omputing the eigenvalues of a matrix.

holdertransformationsandGivensrotations. WeusetheGram-Shmidtpro-

ess. The Gram-Shmidt proess is applied to olumns of the matrix

A

^of

full olumnrank, using the inner produt

h u, v i = u ^∗ v

^.

Denition 2.12. Let

V

^{be an}

n

dimensional vetor spae. A projetion of a vetor

x ∈ V

^{onto the subspae spanned by a vetor}

b

^{is the vetor}

u

intothe same diretionas

b

^{with length}

| u | = | x | cos θ

^{, where}

θ

^{is the angle}

between the vetors

x

^and

b

^{. Beause}

cos θ = x · b

| x || b |

and beause

u

^{is in the diretionof}

b

^{we get}

u = | x | x · b

| x || b | b

| b | .

Hene wean dene: a projetion of

a

^{ontothe subspae spanned by}

e

^is

proj

e a = h e, a i h e , e i e,

where the inner produt

h·i

^{is dened as}

x ^∗ x

^.

x y

y−projection

vector v

x−projection of v

of v

Orthonormalize the olumns of

A = [a ₁ , a ₂ . . . a _n ]

^.

u ₁

⁼

a ₁

^,

e ₁ = _|| ^u _u ¹

1 ||

^,

u ₂

⁼

a ₂ −

^proj

e1 a ₂

^,

e ₂ = _|| ^u u ² ₂ ||

^,

u _n

⁼

a ₃ −

^proj

e1 a ₃ −

^proj

e2 a ₃

^,

e ₃ = _|| ^u u ³ ₃ ||

^,

.

.

.

u _n

⁼

a _n − P n − 1

j=1

^proj

ej a _n

^,

e _n = _|| ^u u ⁿ _n ||

^.

Byrearrangingtheaboveequationssothat the

a _i

^{'sareonthelefthandside,}

and using the fatthat

e _i

^{'sare unit vetors youget:}

a ₁

⁼

h e ₁ , a ₁ i e ₁

^,

a ₂

⁼

h e ₁ , a ₂ i e ₁ + h e ₂ , a ₁ i e ₂

^,

a ₃

⁼

h e ₁ , a ₃ i e ₁ + h e ₂ , a ₃ i e ₂ + h e ₃ , a ₃ i e ₃

^,

.

.

.

a _n

⁼

P n

j=1 h e _j , a _n i e _j

^.

This an bewritten inthe matrix form

A = QR,

where

Q = [e 1 e ₂ . . . e _n ]

^and

R =











h e ₁ , a ₁ i h e ₁ , a ₂ i h e ₁ , a ₃ i . . . 0 h e ₂ , a ₂ i h e ₂ , a ₃ i 0 0 h e ₃ , a ₃ i . . .

.

.

.

.

.

. .

.

.









 .

Example 2.13. Compute aQR-deompositionfor the matrix

A

^,when

A =





2 1 3

− 1 0 7 0 − 1 − 1



 .

The olumns of A are:

a ₁ =



 2

− 1 0



 , a ₂ =



 1 0

− 1



 , a ₃ =



 3 7

− 1



 .

Then use the Gram-Shmidtproess toorthonormalizethe vetors:

q ₁ = a ₁

|| a ₁ || =





√ 2 5

− ^{√ 1} ₅ 0



 ,

q ₂ =

A ₂ − h a ₂ , q ₁ i h q ₁ , q ₁ i

1

|| a ₂ || =







√ 1 2 30

√ 30

− ^{√ 5} ₃₀





 ,

q ₃ =

a ₃ − h a ₃ , q ₁ i

h q ₁ , q ₁ i − h a ₃ , q ₂ i h q ₂ , q ₂ i

1

|| a ₃ || =







√ 1 2 6

√ 6

√ 1 6





 .

Thus you obtainthe orthogonal matrix

Q

^:

A =







√ 2 5

√ 1 30

√ 1 6

− ^{√ 1} ₅ ^{√ 2} ₃₀ ^{√ 2} ₆ 0 − ^{√ 5} ₃₀ ^{√ 1} ₆





 .

The matrix

R

^is

R =





h a ₁ , q ₁ i h a ₂ , q ₁ i h a ₃ , q ₁ i 0 h a ₂ , q ₂ i h a ₃ , q ₃ i 0 0 h a ₃ , q ₃ i



 =







√ 5 ^{√ 2} ₅ − ^{√ 1} ₅ 0 ^{√ 6} ₃₀ ^{√ 22} ₃₀ 0 0 ^{√ 16} ₆





 .

The QR-deompositionis:

A = QR =





2 1 3

− 1 0 7 0 − 1 − 1



 =







√ 2 5

√ 1 30

√ 1 6

− ^{√ 1} ₅ ^{√ 2} ₃₀ ^{√ 2} ₆ 0 − ^{√ 5} ₃₀ ^{√ 1} ₆













√ 5 ^{√ 2} ₅ − ^{√ 1} ₅ 0 ^{√ 6} ₃₀ ^{√ 22} ₃₀ 0 0 ^{√ 16} ₆





 .

TheQR-deompositionanbeomputedmoregenerallyforan

m × n

^matrix

A

^{,as long as}

m ≤ n

^.

QR-deomposition in MATLAB

QR-deomposition is oered as a MATLAB funtion qr. Here is a brief ex-

ampleonhowtosolvelinearsystemsusingQR-deomposition. Thefuntion

triusolve is a self-made funtion todo the bak substitution; reating one

is anexerise task.

Example 2.14. Here is an example on how to use QR-deomposition to

solvea system of linear equations.

A =

3 -5 7

0 4 5

-6 -9 -8

>>x = ones(3,1);

>>b = A*x

b =

5

9

-23

>>[q r℄ = qr(A)

>> ba = q'*b;

>> xs = triusolve(r,ba);

>>norm(A*xs-b)

ans =

6.6465e-15

2.4.4 Singular value deomposition

Every

m × n

^{matrix with omplex entries an bepresented as the produt}

A = USV ^∗ ,

where

U ∈ C ^{m × m}

and

V ∈ C ^{n × n}

are orthogonal matries, and

S ∈ C ^{m × n}

is a diagonal matrix with entries sorted by magnitude. If the matrix

A

^is

invertible, the inverse is

A ^{− 1} = VS ^{− 1} U ^∗ .

Thisiseasytoompute, astheinverse ofadiagonalmatrix isjustadiagonal

matrix with inverses of the originaldiagonalelements.

MATLAB uses the QR-algorithm to ompute the singular value deompo-

sition. If one wishes to ompute it manually, the following proedure is

propably the easiest:

1. Find the eigenvalues and orthonormalized eigenvetors of

A ^∗ A

^,i.e.,

A ^∗ A = VΛV ^∗ .