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)

 

sqrt(5), sqrt(-5)

 

ans = 2.2361

ans = 0 + 2.2361i

 

% A quadratic equation, a*x^2 + b*x + c = 0, has two solutions (roots):

% 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

 

c1 = round(1.5), c2 = round(1.4999999), c3 = round(-1.5), c4 = round(-1.4999999)

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 = rand(1), a2 = rand(1), a3 = rand(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:

 

         function m-files:

 

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

 

 

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'.

 

 

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

 

 

% 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