# Calcolo numerico - il Tutorial Matlab

Anteprima

### ESTRATTO DOCUMENTO

2

Function Description

Inverse cosine

acos Control axis scaling and appearance

axis Create character array

char Cholesky factorization

chol Cosine function

cos Vector cross product

cross Determinant

det Diagonal matrices and diagonals of a matrix

diag Convert to double precision

double Eigenvalues and eigenvectors

eig Identity matrix

eye Filled 2-D polygons

fill Round towards zero

fix Flip matrix in left/right direction

fliplr Floating point operation count

flops Grid lines

grid Hadamard matrix

hadamard Hilbert matrix

hilb Hold current graph

hold Matrix inverse

inv True for empty matrix

isempty Graph legend

legend Length of vector

length Linearly spaced vector

linspace Convert numerical values to logical

logical Magic square

magic Largest component

max Smallest component

min Matrix or vector norm

norm Null space

null Convert numeric array into cell array

num2cell Convert number to string

num2str Ones array

ones Pascal matrix

pascal Linear plot

plot Convert roots to polynomial

poly Evaluate polynomial

polyval Uniformly distributed random numbers

rand Normally distributed random numbers

randn Matrix rank

rank Reduced row echelon form

reff Remainder after division

rem Change size

reshape Find polynomial roots

roots Sine function

sin Size of matrix

size Sort in ascending order

sort 3

Symbolic substitution

subs Construct symbolic bumbers and variables

sym Start a stopwatch timer

tic Graph title

title Read the stopwatch timer

toc Tioeplitz matrix

toeplitz Extract lower triangular part

tril Extract upper triangular part

triu Vandermonde matrix

vander Variable length input argument list

varargin Zeros array

zeros

The purpose of this section is to demonstrate how to create and transform vectors and matrices in

## MATLAB.

This command creates a row vector

a = [1 2 3]

a = 1 2 3

Column vectors are inputted in a similar way, however, semicolons must separate the components

of a vector

b = [1;2;3]

b = 1

2

3

The quote operator ' is used to create the conjugate transpose of a vector (matrix) while the dot-

quote operator .' creates the transpose vector (matrix). To illustrate this let us form a complex

vector a + i*b' and next apply these operations to the resulting vector to obtain

(a+i*b')'

ans =

1.0000 - 1.0000i

2.0000 - 2.0000i

3.0000 - 3.0000i

while 4

(a+i*b').'

ans =

1.0000 + 1.0000i

2.0000 + 2.0000i

3.0000 + 3.0000i

Command length returns the number of components of a vector

length(a)

ans = 3

The dot operator. plays a specific role in MATLAB. It is used for the componentwise application

of the operator that follows the dot operator

a.*a

ans = 1 4 9

The same result is obtained by applying the power operator ^ to the vector a

a.^2

ans = 1 4 9

Componentwise division of vectors a and b can be accomplished by using the backslash operator

\ together with the dot operator .

a.\b'

ans = 1 1 1

For the purpose of the next example let us change vector a to the column vector

a = a'

a = 1

2

3

The dot product and the outer product of vectors a and b are calculated as follows

dotprod = a'*b 5

dotprod =

14

outprod = a*b'

outprod =

1 2 3

2 4 6

3 6 9

The cross product of two three-dimensional vectors is calculated using command cross. Let the

vector a be the same as above and let

b = [-2 1 2];

Note that the semicolon after a command avoids display of the result. The cross product of a and

b is

cp = cross(a,b)

cp = 1 -8 5

The cross product vector cp is perpendicular to both a and b

[cp*a cp*b']

ans = 0 0

We will now deal with operations on matrices. Addition, subtraction, and scalar multiplication are

defined in the same way as for the vectors.

This creates a 3-by-3 matrix

A = [1 2 3;4 5 6;7 8 10]

A = 1 2 3

4 5 6

7 8 10

Note that the semicolon operator ; separates the rows. To extract a submatrix B consisting of

rows 1 and 3 and columns 1 and 2 of the matrix A do the following

B = A([1 3], [1 2])

B = 1 2

7 8

To interchange rows 1 and 3 of A use the vector of row indices together with the colon operator

C = A([3 2 1],:) 6

C = 7 8 10

4 5 6

1 2 3

The colon operator : stands for all columns or all rows. For the matrix A from the last example

the following command

## A(:)

ans = 1

4

7

2

5

8

3

6

10

creates a vector version of the matrix A. We will use this operator on several occasions.

To delete a row (column) use the empty vector operator [ ]

A(:, 2) = []

A = 1 3

4 6

7 10

Second column of the matrix A is now deleted. To insert a row (column) we use the technique for

creating matrices and vectors

A = [A(:,1) [2 5 8]' A(:,2)]

A = 1 2 3

4 5 6

7 8 10

Matrix A is now restored to its original form.

Using MATLAB commands one can easily extract those entries of a matrix that satisfy an impsed

condition. Suppose that one wants to extract all entries of that are greater than one. First, we

define a new matrix A

A = [-1 2 3;0 5 1]

A = -1 2 3

0 5 1 7

Command A > 1 creates a matrix of zeros and ones

A > 1

ans = 0 1 1

0 1 0

with ones on these positions where the entries of A satisfy the imposed condition and zeros

everywhere else. This illustrates logical addressing in MATLAB. To extract those entries of the

matrix A that are greater than one we execute the following command

A(A > 1)

ans = 2

5

3

The dot operator . works for matrices too. Let now

A = [1 2 3; 3 2 1] ;

The following command

## A.*A

ans = 1 4 9

9 4 1

computes the entry-by-entry product of A with A. However, the following command

## A*A

¨??? Error using ==> *

Inner matrix dimensions must agree.

generates an error message.

Function diag will be used on several occasions. This creates a diagonal matrix with the diagonal

entries stored in the vector d

d = [1 2 3];

D = diag(d)

D = 1 0 0

0 2 0

0 0 3 8

To extract the main diagonal of the matrix D we use function diag again to obtain

d = diag(D)

d = 1

2

3

What is the result of executing of the following command?

diag(diag(d));

In some problems that arise in linear algebra one needs to calculate a linear combination of

several matrices of the same dimension. In order to obtain the desired combination both the

coefficients and the matrices must be stored in cells. In MATLAB a cell is inputted using curly

braces{ }. This

c = {1,-2,3}

c = [1] [-2] [3]

is an example of the cell. Function lincomb will be used later on in this tutorial.

function M = lincomb(v,A)

% Linear combination M of several matrices of the same size.

% Coefficients v = {v1,v2,…,vm} of the linear combination and the

% matrices A = {A1,A2,...,Am} must be inputted as cells.

m = length(v);

[k, l] = size(A{1});

M = zeros(k, l);

for i = 1:m

M = M + v{i}*A{i};

end

! "

MATLAB has several tool needed for computing a solution of the system of linear equations.

Let A be an m-by-n matrix and let b be an m-dimensional (column) vector. To solve

the linear system Ax = b one can use the backslash operator \ , which is also called the left

division. 9

1. Case m = n

In this case MATLAB calculates the exact solution (modulo the roundoff errors) to the system in

question.

Let

A = [1 2 3;4 5 6;7 8 10]

A = 1 2 3

4 5 6

7 8 10

and let

b = ones(3,1);

Then

x = A\b

x = -1.0000

1.0000

0.0000

In order to verify correctness of the computed solution let us compute the residual vector r

r = b - A*x

r =

1.0e-015 *

0.1110

0.6661

0.2220

Entries of the computed residual r theoretically should all be equal to zero. This example

illustrates an effect of the roundoff erros on the computed solution

.

2. Case m > n

If m > n, then the system Ax = b is overdetermined and in most cases system is inconsistent. A

solution to the system Ax = b, obtained with the aid of the backslash operator \ , is the least-

squares solution.

Let now

A = [2 –1; 1 10; 1 2];

and let the vector of the right-hand sides will be the same as the one in the last example. Then 10

x = A\b

x = 0.5849

0.0491

The residual r of the computed solution is equal to

r = b - A*x

r = -0.1208

-0.0755

0.3170

Theoretically the residual r is orthogonal to the column space of A. We have

r'*A

ans =

1.0e-014 *

0.1110

0.6994

3. Case m < n

If the number of unknowns exceeds the number of equations, then the linear system is

underdetermined. In this case MATLAB computes a particular solution provided the system is

consistent. Let now

A = [1 2 3; 4 5 6];

b = ones(2,1);

Then

x = A\b

x = -0.5000

0

0.5000

A general solution to the given system is obtained by forming a linear combination of x with the

columns of the null space of A. The latter is computed using MATLAB function null

z = null(A)

z = 0.4082

-0.8165

0.4082 11

Suppose that one wants to compute a solution being a linear combination of x and z, with

coefficients 1 and –1. Using function lincomb we obtain

w = lincomb({1,-1},{x,z})

w = -0.9082

0.8165

0.0918

The residual r is calculated in a usual way

r = b - A*w

r =

1.0e-015 *

-0.4441

0.1110

# $

The built-in function rref allows a user to solve several problems of linear algebra. In this section

we shall employ this function to compute a solution to the system of linear equations and also to

find the rank of a matrix. Other applications are discussed in the subsequent sections of this

tutorial.

Function rref takes a matrix and returns the reduced row echelon form of its argument. Syntax of

the rref command is or

B = rref(A) [B, pivot] = rref(A)

The second output parameter pivot holds the indices of the pivot columns.

Let

A = magic(3); b = ones(3,1);

A solution x to the linear system Ax = b is obtained in two steps. First the augmented matrix of

the system is transformed to the reduced echelon form and next its last column is extracted

[x, pivot] = rref([A b])

x = 1.0000 0 0 0.0667

0 1.0000 0 0.0667

0 0 1.0000 0.0667

pivot =

1 2 3 12

x = x(:,4)

x = 0.0667

0.0667

0.0667

The residual of the computed solution is

b - A*x

ans = 0

0

0

Information stored in the output parameter pivot can be used to compute the rank of the matrix A

length(pivot)

ans = 3

% &

MATLAB function inv is used to compute the inverse matrix.

Let the matrix A be defined as follows

A = [1 2 3;4 5 6;7 8 10]

A = 1 2 3

4 5 6

7 8 10

Then

B = inv(A)

B = -0.6667 -1.3333 1.0000

-0.6667 3.6667 -2.0000

1.0000 -2.0000 1.0000

In order to verify that B is the inverse matrix of A it sufficies to show that A*B = I and

B*A = I, where I is the 3-by-3 identity matrix. We have 13

## A*B

ans =

1.0000 0 -0.0000

0 1.0000 0

0 0 1.0000

In a similar way one can check that B*A = I.

The Pascal matrix, named in MATLAB pascal, has several interesting properties. Let

A = pascal(3)

A = 1 1 1

1 2 3

1 3 6

Its inverse B

B = inv(A)

B = 3 -3 1

-3 5 -2

1 -2 1

is the matrix of integers. The Cholesky triangle of the matrix A is

S = chol(A)

S = 1 1 1

0 1 2

0 0 1

Note that the upper triangular part of S holds the binomial coefficients. One can verify easily that

A = S'*S.

Function rref can also be used to compute the inverse matrix. Let A is the same as above. We

create first the augmented matrix B with A being followed by the identity matrix of the same size

as A. Running function rref on the augmented matrix and next extracting columns four through

six of the resulting matrix, we obtain

B = rref([A eye(size(A))]);

B = B(:, 4:6)

B = 3 -3 1

-3 5 -2

1 -2 1 14

To verify this result, we compute first the product A *B

## A*B

ans = 1 0 0

0 1 0

0 0 1

and next B*A

## B*A

ans = 1 0 0

0 1 0

0 0 1

This shows that B is indeed the inverse matrix of A.

' (

In some applications of linear algebra knowledge of the determinant of a matrix is required.

MATLAB built-in function det is designed for computing determinants.

Let

A = magic(3);

Determinant of A is equal to

det(A)

ans =

-360

One of the classical methods for computing determinants utilizes a cofactor expansion. For more

details, see e.g., [2], pp. 103-114. entry of the matrix A

Function ckl = cofact(A, k, l) computes the cofactor ckl of the a kl

function ckl = cofact(A,k,l)

% Cofactor ckl of the a_kl entry of the matrix A.

[m,n] = size(A);

if m ~= n

error('Matrix must be square') 15

end

B = A([1:k-1,k+1:n],[1:l-1,l+1:n]);

ckl = (-1)^(k+l)*det(B);

Function d = mydet(A) implements the method of cofactor expansion for computing

determinants

function d = mydet(A)

% Determinant d of the matrix A. Function cofact must be

% in MATLAB's search path.

[m,n] = size(A);

if m ~= n

error('Matrix must be square')

end

a = A(1,:);

c = [];

for l=1:n

c1l = cofact(A,1,l);

c = [c;c1l];

end

d = a*c;

Let us note that function mydet uses the cofactor expansion along the row 1 of the matrix A.

Method of cofactors has a high computational complexity. Therefore it is not recommended for

computations with large matrices. Its is included here for pedagogical reasons only. To measure a

computational complexity of two functions det and mydet we will use MATLAB built-in

function flops. It counts the number of floating-point operations (additions, subtractions,

multiplications and divisions). Let

A = rand(25);

be a 25-by-25 matrix of uniformly distributed random numbers in the interval ( 0, 1 ). Using

function det we obtain

flops(0)

det(A)

ans =

-0.1867

flops

ans = 10100

For comparison, a number of flops used by function mydet is

flops(0) 16

mydet(A)

ans =

-0.1867

flops

ans = 223350

The adjoint matrix adj(A) of the matrix A is also of interest in linear algebra (see, e.g., [2],

p.108).

function B = adj(A)

% Adjoint matrix B of the square matrix A.

[m,n] = size(A);

if m ~= n

error('Matrix must be square')

end

B = [];

for k = 1:n

for l=1:n

B = [B;cofact(A,k,l)];

end

end

B = reshape(B,n,n);

The adjoint matrix and the inverse matrix satisfy the equation

-1

A = adj(A)/det(A)

(see [2], p.110 ). Due to the high computational complexity this formula is not recommended for

computing the inverse matrix.

)

The 2-norm (Euclidean norm) of a vector is computed in MATLAB using function norm.

Let

a = -2:2

a = -2 -1 0 1 2

The 2-norm of a is equal to

twon = norm(a) 17

twon =

3.1623

With each nonzero vector one can associate a unit vector that is parallel to the given vector. For

instance, for the vector a in the last example its unit vector is

unitv = a /twon

unitv =

-0.6325 -0.3162 0 0.3162 0.6325

θ

The angle between two vectors a and b of the same dimension is computed using the formula

= arccos(a.b/||a|| ||b||),

where a.b stands for the dot product of a and b, ||a|| is the norm of the vector a and arccos is the

inverse cosine function.

Let the vector a be the same as defined above and let

b = (1:5)'

b = 1

2

3

4

5

Then

angle = acos((a*b)/(norm(a)*norm(b)))

angle =

1.1303

Concept of the cross product can be generalized easily to the set consisting of n -1 vectors in the

n

n-dimensional Euclidean space . Function crossprod provides a generalization of the

MATLAB function cross.

function cp = crossprod(A)

% Cross product cp of a set of vectors that are stored in columns of A.

[n, m] = size(A);

if n ~= m+1

error('Number of columns of A must be one less than the number of

rows')

I contenuti di questa pagina costituiscono rielaborazioni personali del Publisher cecilialll di informazioni apprese con la frequenza delle lezioni di Calcolo numerico e studio autonomo di eventuali libri di riferimento in preparazione dell'esame finale o della tesi. Non devono intendersi come materiale ufficiale dell'università Napoli Federico II - Unina o del prof D'Amore Luisa.

Acquista con carta o conto PayPal

Scarica il file tutte le volte che vuoi

Paga con un conto PayPal per usufruire della garanzia Soddisfatto o rimborsato