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