%------------------------------------------------------------------------
%%%               ENGINEERING MATHEMATICS III - MATH 2Z03
%%%                               LAB #1
%%%                       INTRODUCTION TO MATLAB
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Covers:
% - "Numerical Mathematics" by M. Grasselli and D. Pelinovsky,
%    Sections 1.3 - 1.6
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Instructions:
% - Execute the commands given below in MATLAB. You can do this in
%   one of the following three ways:
%   (i) retype the commands one by one at the MATLAB command line (>>),
%   (ii) copy and paste the commands one by one to the MATLAB
%        command line (>>),
%   (iii) run all the commands at once by calling the present file
%         as a script, i.e. >> m2z03lab
% - Analyze the output and compare it with the comments embedded in
%   the file below.
%------------------------------------------------------------------------

clc; close all; clear all;

%------------------------------------------------------------------------
% Section 1.3: Matrix Operations
disp(' -----------------------------------------------------------')
disp(' Section 1.3: Matrix Operations ')

% matrix addition and subtraction 
% note consistency of the dimensions in all examples below;
A=[2 -3 0; 1 2 3; 0 0 1];
B=[1 0 0; 13 12 11; -1 -1 3];
C=[5 -5; 2 1]; D=[1 15; -3 4];
A+B

C-D

% matrix multiplication (note consistency of the dimensions);
A=[5 2 3; 1 -2 6];
B=[0 1 -1 3; 5 1 0 2; 4 -3 2 -1];
C=A*B

A=[1 0 -1;5 -3 6; 7 3 1];
B=[2 2 1; -1 4 1; 0 0 2];
A*B

B*A


% matrix-vector multiplication;
A=[5 2 3; 1 -2 6];
x=[-1 1 0]';
b=A*x

y=[2 -17];
c=y*A


% the dot (scalar, inner) product;
a=[7;2;-9]; b=[0;-1;3];
a'*b

dot(a,b)

% the outer (tensor) product;
a*b'

% the cross product;
a=[1;3;-4]; b=[2;-1;0]; c =[5;7;1];
[cross(a,b),cross(b,a)]

% checking some identities of the vector calculus;
[cross(a,cross(b,c)),cross(cross(a,b),c)]

[dot(a,cross(a,b)),dot(b,cross(a,b))]

% matrix exponentiation (required matrix inverses for negative
% exponents);
A=[1 0; -1 2];
A^3

A^(-2)

% componentwise operation (note that algebraic operators are
% preceded with a dot);
A=[1 0; -1 2]; B=[-1 3; 4 1];
A.*B

A./B

A.^2

B.^A

disp('  Press any key to continue ... ')
pause
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Section 1.4: Built-in Functions
disp(' -----------------------------------------------------------')
disp(' Section 1.4: Built-in Functions ')

% Use command 'help' to obtain information about any MATLAB function
% e.g. componentwise calculation of square root
help sqrt

% Compare:
A = [1 4 -9]
sqrt(A)
A.^(1/2)

% or matrix square root
help sqrtm

%%% Finding matrix inverse
help inv

% Non-singular case
A = [1 3; -1 2]
inv(A)

% Singular case (with a warning message)
B = [-1 3; 2 -6]
inv(B)

%%% Matrix division
% Using matrix inverse command
A = [1 3; -1 2]
B = [-1 4; 5 -6]
X = A * inv(B)

% Using matrix division
X = A / B

disp('  Press any key to continue ... ')
pause
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Section 1.5: Programming with MATLAB
disp(' -----------------------------------------------------------')
disp(' Section 1.5: Programming with MATLAB ')

% Function marks.m
%------------------------------------------------------------------------
%function [m,s] = marks(x)
% % [m,s] = marks(x)
% % x is a vector, m and s are scalars
% % the function computes the mean value m and standard deviation s
% % of the data sample x

%n = length(x);
%m = sum(x)/n; % built-in function for the sum of coponents of a vector
%s = sqrt(sum((x-m).^2)/n); % componentwise operations
%------------------------------------------------------------------------

% Type it (or copy and paste) as a new M-file removing the sign '%'
% at the beginning of every line; Save it as marks.m

% Input data for the function marks.m
x = [80 70 85 65 70 75 90];

% Calling the function marks.m
[m,s] = marks(x)


%%% Branch statements (if .. else .. elseif .. end)
% This script calculates final letter grade based on marks enterd by user

% Invitation for the user to create a vector x as input data
x = input('Enter your term marks as a vector: ');

% Calling the function marks.m
[m,s] = marks(x);

% Presenting the results returned by the function marks.m 
disp('  Your numerical average is: ')
disp(m)

% Branch statements depending on the value of m
if m >= 80 % condition
    disp('Your final grade is A - well done!') % message to show if condition holds
elseif (70 <= m)&&(m < 80)
    disp('Your final grade is B')
elseif (60 <= m)&&(m < 70)
    disp('Your final grade is C')
elseif (50 <= m)&&(m < 60)
    disp('Your final grade is D')
else
    disp('Sorry, but you failed the course!')
end

% Enter any vector of marks, e.g. [50 61 70 56 80 75]


%%% Loops (for .. end)
% This script calculates the sum of squares for the first N integers

% Number of the first N integers
N = 1000

% Initial value of the sum
s = 0;

% Loop
for x = 1:N % vectorization of x
    s = s + x^2; % summing N times
end

% Result output
s

% Compare with usage of the colon operator ':'
sum([1:1000].^2)


%%% Loops (while .. end)
% This script computes how many terms k of a geometric progression starting
% at a0 with ratio q need to be added until you reach a number grater than
% a certain threshold M

% Initial data
a0 = 1
q = 2
M = 40 * 1000 * 5 * 400 * 10^9

% First term value and its number
x(1) = a0;
k = 1;

% Loop
while sum(x) < M % condition to check (loop ends when it is faulse)
    x(k+1) = q * x(k);
    k = k + 1;
end

% Result output
k

disp('  Press any key to continue ... ')
pause
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Section 1.6: Graphics and Data Files
disp(' -----------------------------------------------------------')
disp(' Section 1.6: Graphics and Data Files ')

%%% Command 'plot'
% Creates a two-dimensional plot of one vector versus another

% Vectorisation using colon operator (argument)
x = 0:0.2:6

% Creating another vector (function) on the basis of x vectorisation
y = exp(-x).*sin(pi*x)

% See result on Figure 1
plot(x,y);

disp('  Press any key to continue ... ')
pause

% To improve the resolution we will refine resolution in the first
% line. Compare the results in Figures 1 and 2
x = 0:0.02:6;
y = exp(-x).*sin(pi*x);
figure; % opening new window for new figure
plot(x,y);

disp('  Press any key to continue ... ')
pause

% Plotting two-dimensional curves given in the parametric form
t = 0:0.01:2*pi; % parameter t
a = 4; b = 2 * sqrt(3); % coefficients
x = a * cos(t) + 2; % x(t)
y = b * sin(t); % y(t)
figure; plot(x,y); % see Figure 3
grid on; % showing grid on graph

disp('  Press any key to continue ... ')
pause

% Embelishments: different line colors, styles, legends
x = 0:0.01:5; % x vectorization
y = sin(x); % y(x)
figure; plot(x,y,'--r'); % Figure 4: red dashed line

z = cos(x); % z(x)
hold on; % command that allows one to place new graphs in the same figure
plot(x,z); % Figure 4: regular style (blue solid line)

legend('sin(x)', 'cos(x)'); % adding legend to current figure
title('SIN(X) and COS(X)'); % adding title to current figure

disp('  Press any key to continue ... ')
pause

%%% Command 'polar'
% Creates a two-dimensional plot of a curve in polar coordinates

t = 0:0.01:2*pi; % parameter t
a = 4; b = 2 * sqrt(3); % coefficients
r = (b^2/a)./(1-sqrt(1-b^2/a^2)*cos(t)); % r(t)
figure; polar(t,r); % see Figure 5

disp('  Press any key to continue ... ')
pause

%%% Command 'plot3'
% Creates three-dimensional curves

z = 0:0.1:10; % z vectorization
x = exp(-0.2*z).*cos(pi*z); % x(z)
y = exp(-0.2*z).*sin(pi*z); % x(z)
figure; plot3(x,y,z); % see Figure 6

disp('  Press any key to continue ... ')
pause

%%% Command 'surf'
% Generates surfaces of the form z = f(x,y)

x = -2:0.1:2; % x vectorization
y = -1:0.1:1; % y vectorization
[X,Y] = meshgrid(x,y); % creating mesh
z = X.*Y.*exp(-X.^2-Y.^2); % z = f(x,y)
figure; surf(x,y,z); % see Figure 7

disp('  This is the end of the lab tutorial #1. ')