CHAPTER 1: INTRODUCTION INTO MATLAB COMPUTING

Lecture 1.3: Mathematical functions and M-files

Elementary mathematical functions:

·         sqrt: square root function (MATLAB operates both with real and complex numbers universally)

ans =    2.2361

ans =    0 + 2.2361i

% x = (-b +/- sqrt(b^2 – 4*a*c))/(2*a)

a = 1; b = 1; c = 10;

x1 = (-b + sqrt(b^2-4*a*c))/(2*a)

x2 = (-b - sqrt(b^2-4*a*c))/(2*a)

% the two solutions are complex since b^2 < 4*a*c

a = 1; b = 10; c = 1;

x1 = (-b + sqrt(b^2-4*a*c))/(2*a)

x2 = (-b - sqrt(b^2-4*a*c))/(2*a)

% the two roots are real since b^2 > 4*a*c

x1 =  -0.5000 + 3.1225i

x2 =  -0.5000 - 3.1225i

x1 =   -0.1010

x2 =   -9.8990

·         sin, cos, tan: trigonometric functions

·         asin, acos, atan: inverse trigonometric functions

·         sinh, cosh, tanh: hyperbolic functions

·         asinh, acosh, atanh: inverse hyperbolic functions

·         exp, log, log10: exponential functions

a1 = cos(3-4*i) % MATLAB functions work for complex arguments

a2 = exp(7+5*i)

a3 = tan(0:1:5) % MATLAB functions work for vector arguments

a4 = cosh([ 0,1; -1,-2; 0.1,0.2])

% MATLAB functions work for matrix arguments, such that

% the output of MATLAB functions reproduces the input structure

% Note that the matrix-vector calls to MATLAB functions are more efficient

% computationally than pointwise calls for each element of an array

a1 =

-27.0349 + 3.8512i

a2 =

3.1107e+002 -1.0516e+003i

a3 =

0    1.5574   -2.1850   -0.1425    1.1578   -3.3805

a4 =

1.0000    1.5431

1.5431    3.7622

1.0050    1.0201

b1 = acos(0.5) % the argument of the standard "acos" is between [-1,1]

b2 = acos(5) % the MATLAB "acos" is analytically continued beyond [-1,1]

b3 = acosh(5) % imaginary part of acos(5) is the same as real part of acosh(5)

b4 = acosh(0.5) % vice verse

b5 = log(exp(5)) % log(x) is the logarithm of base e

b6 = log10(1000) % log10(x) is the logarithm of base 10

format long; b7 = atan(inf) % the range of atan(x) is between [-pi/2,pi/2]

pi/2, format short;

b1 =    1.0472

b2 =    0 - 2.2924i

b3 =    2.2924

b4 =    0 + 1.0472i

b5 =    5

b6 =    3.0000

b7 =   1.57079632679490

ans =  1.57079632679490

·         real, imag, conj: real, imaginary parts of a complex number, and the complex conjugate number

·         abs, angle:  absolute value of a complex number, and the phase angle of a complex number

z = 1 + 2*i; % a complex number

r = abs(z), theta = angle(z) % polar form of a complex number

comp1 = real(z) - r*cos(theta) % Euler formulas for complex exponentials

comp2 = imag(z) - r*sin(theta)

r =    2.2361

theta =    1.1071

comp1 = -2.2204e-016

comp2 =     0

·         round: rounding to the nearest integer

·         floor: rounding to the smallest integer

·         ceil: rounding to the largest integer

·         fix: rounding to the smallest integer if positive and to the largest integer if negative

·         sign: 1 if positive, 0 if zero, and –1, if negative

d1 = floor(1.9999), d2 = floor(-1.0001), d3 = ceil(1.0001), d4 = ceil(-1.99999)

e1 = fix(1.99999), e2 = fix(-1.99999), e3 = sign(1.01), e4 = sign(-1.01)

c1 =     2

c2 =     1

c3 =    -2

c4 =    -1

d1 =     1

d2 =    -2

d3 =     2

d4 =    -1

e1 =     1

e2 =    -1

e3 =     1

e4 =    -1

·         mod,rem: remainders upon division

% mod(x,y) = x – y*ceil(x/y)

% rem(x,y) = x – y*fix(x/y)

% rem(x,y) = mod(x,y) if sign(x) = sign(y)

f1 = mod(5,3), f2 = rem(5,3)

f3 = mod(-5,3), f4 = rem(-5,3)

f5 = mod(5,-3), f6 = rem(5,-3)

f7 = mod(-5,-3), f8 = rem(-5,-3)

f1 =     2

f2 =     2

f3 =     1

f4 =    -2

f5 =    -1

f6 =     2

f7 =    -2

f8 =    -2

·         factorial: factorial of an integer number

·         factor: factorize an integer number into prime numbers

·         isprime: 1 if a number is a prime number and 0 if not

g1 = factorial(5), g2 = factorial(15), g3 = factor(150), g4 = isprime(150)

g1 =   120

g2 =  1.3077e+012

g3 =     2     3     5     5

g4 =     0

g5 = factorial(5.1),

% the functions are not continued for non-integer values

??? Error using ==> factorial

N must be a positive integer

% Remove the numbers divisible by 3 from a vector:

x = [-9, 0, 2, 5, 3, 7, 2, 0, 6, 12, 214342124187]

disp(sprintf(' %12.0f ',x));

m = 0;

for k = 1 : length(x) % looping through the vector x

if mod(x(k),3) ~= 0 % comparison if x(k) is not divisible by 3

m = m + 1;

y(m) = x(k); % saving x(k) into a dynamically created vector y

end

end

x = y

x =

1.0e+011 *

-0.0000         0    0.0000    0.0000    0.0000    0.0000    0.0000         0    0.0000    0.0000    2.1434

-9             0             2             5             3             7             2             0             6            12  214342124187

x =

2     5     7     2

Standard API (Application Programming Interface)_functions:

·         sort: the function reorders elements of a vector to ascending order

x = [2 1 4 10 7 3 1 ]

y = sort(x) % both row-vectors and column-vectors can be sorted

[y,index] = sort(x); % index show how elements of y appear in x: y = x(ind)

x(ind)

x =     2     1     4    10     7     3     1

y =     1     1     2     3     4     7    10

ans =   1     1     2     3     4     7    10

·         sum:  the function summates all elements of a vector

sum(x), sum(y')

A = [ 1,2,3; -4,-5,-6] % a two-dimensional array (matrix)

sum(A) % returns the row-vector for sums of each column of A

ans =    28

ans =    28

A =     1     2     3

-4    -5    -6

ans =

-3    -3    -3

·         max,min: find the maximum, minimum element of a vector

[maxX,indX] = max(x) % maxX - the maximal element, indX - its position in x

[minX,indX] = min(x) % if the minimum element is multiple, indX is the index of the first element

[maxA,indA] = max(A) % returns the row vectors for max of each column of A

maxX =    10

indX =     4

minX =     1

indX =     2

maxA =     1     2     3

indA =     1     1     1

·         rand(n): generate a square matric (n by n) of random numbers, each random number is in [0,1]

a1 =    0.4565

a2 =    0.0185

a3 =    0.8214

·         rand('seed',k): initialization of the random number generator by a seed number k (>=1)

NB: If the seed number is the same, the sequence of random numbers becomes the same. In order to generate a different sequence of random numbers, the seed number is chosen according to a random factor (e.g. by count of the day, time in seconds, etc.)

c = clock % returns a row-vector of length 6

k = c(2)*c(3)*c(4)*c(5)*c(6) % this product has 3*10^7 combinations and changes every second during the year

rand('seed',k); % initialization of the random number generator depends on the time of computations

for n = 1 : 10

m = ceil(12*rand(1)); % randomly choose one the 12 integer numbers between 1 and 12 for 10 times

fprintf('the random number is: %3.0f\n',m);

end

c =  1.0e+003 *

2.0010    0.0120    0.0280    0.0110    0.0190    0.0023

k =  1.5934e+005

the random number is:  11

the random number is:   4

the random number is:  10

the random number is:   6

the random number is:   3

the random number is:   9

the random number is:   1

the random number is:   6

the random number is:   8

the random number is:   5

·         eval('stringName'): execute the command prescribed by the string "stringName"

FunctionName = input('Type a function name: ','s');

% an user is supposed to type a name of valid MATLAB function

x = 0 : 0.2 : 1;

strCommand = [ FunctionName, '(x)' ]; % concatanetion of two strings

y = eval(strCommand) % evaluate the command, e.g. y = sin(x)

Type a function name: sin

y =

0    0.1987    0.3894    0.5646    0.7174    0.8415

Type a function name: exp

y =

1.0000    1.2214    1.4918    1.8221    2.2255    2.7183

% Example of authomatical creating names of files:

for k = 1 : 1000

y = exp((k-1)*0.1); % the data y is dynamically created for each value of k

if ( k <= 10 )

s = sprintf('save file00%1d y -ascii',k-1); % saving the data y into file "file00k" for k < 10

elseif ( k <= 100 )

s = sprintf('save file0%1d y -ascii',k-1); % saving the data y into file "file0k" for 10 <= k < 100

else

s = sprintf('save file%1d y -ascii',k-1); % saving the data y into file "filek" for 100 <= k < 1000

end

eval(s); % executing the string with the command to save the data, 1000 files are created in the loop

end

M-file scripts and functions:

M-files are useful for long commands and procedures, they can be saved to disk and corrected as many times as needed.

·         script m-files:

• a script m-file corresponds to a main program in other programming languages
• a script m-file can include anything the user would write in the Command Window
• all comments in the beginning of a script can be viewed by on-line help
• echo on: all statements are printed out on the screen while they are being executed
• a script m-file can be called inside other script m-files, all variables in the current working space are shared between the parent and the children m-file

·         function m-files:

• function m-files correspond to subroutines or functions in other programming languages
• each function is saved as a separate m-file
• a function m-file starts with the line: function [y1,y2,…,yN] = functionName(x1,x2,…,xM)
• y1,y2,…,yN – N output values
• x1,x2,…,xM – M input arguments
• functionName – the same as the name of the m-file
• all comments after the first line (declaration of the function) can be viewed by on-line help as a description of the MATLAB function
• a core of the function m-file is a script that uses x1,x2,…,xM, perform any intermediate computations and finally compute y1,y2,…,yN

function [mean,stdev] = mean_stdev(x)

% [mean,stdev] = mean_stdev(x)

% x is a vector, mean and stdev are scalars

% the function computes the mean value and the standard deviation

% of a data sample x

n = length(x); % x is a vector

if ( n ~= 0 ) % if x is not empty vector

mean = sum(x)/n; % compute the mean value of x

varComp = sum((x-mean).^2)/n; % compute the variance of x

stdev = sqrt(varComp); % compute the standard deviation of x

end

• a function m-file can be called in another function m-file, or script m-file, or directly in Command Window, the function shares only variables given as input arguments

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

[mu,st] = mean_stdev(x)

x = 3;

[mu,st] = mean_stdev(x)

mu =    5.5000

st =    2.8723

mu =     3

st =     0

mean_stdev(3)

varComp % the variable "varComp" is defined in the function "mean_stdev",

% but this variable is not accessible in the current working space

ans =

3

??? Undefined function or variable 'varComp'.

• MATLAB function m-file should be generalized for vector and matrix arguments

function y = powerfunction(x,n)

y = x.^n; % the pointwise power ".^" works both for scalars and vectors

x = 3; powerfunction(x,2)

x = [ 1 2 3 4 5 6 ]; powerfunction(x,3)

ans =

9

ans =

1     8    27    64   125   216

• feval: evaluate a function m-file, e.g. [y1,y2,…,yN] = feval(f_name,x1,x2,…,xM), where "f_name" is the string for the function name; x1,x2,…,xM are values of the input arguments, and y1,y2,…,yN are the output values.

% compute a weighted average of a function y = f(x) at three points a,b,c by using the formula:

% f_av = (f(a) + 2 f(b) + f(c))/4

function f_av  = average(f,a,b,c)

a1 = feval(f,a);

b1 = feval(f,b);

c1 = feval(f,c);

f_av = 0.25*(a1 + 2*b1 + c1);

w1 = average('sqrt',1,2,3)

% call to the standard MATLAB function y = sqrt(x)

w2 = average('sin',1,2,3)

% call to the standard MATLAB function y = sin(x)

w1 =

1.3901

w2 =

0.7003