echo on
%%%%% lab2_matlab.m %%%%%

disp(sprintf('**** Exercises for practicing MatLab ****\n'));

disp(sprintf('**** VECTOR NORMS *'));

v=[-1.6 1.2]';	% Vector norms:
%norm(v,p)	% Returns sum(abs(Vi).^p)^(1/p), for i=1..n and any p.
norm(v,1)	% Sum of absolute values of elements of vector V.
norm(v,2)	% Also norm(V). 
norm(v,inf)	% Returns max(abs(Vi)).

disp(sprintf('**** MATRIX NORMS *****'));

A=[2 -1 1; 1 0 1; 3 -1 4];	% Matrix norms:
norm(A,1)	% The 1-norm is largest absolute column sum of A.
norm(A,2)	% The largest singular value of A (same as norm(A)).
norm(A,inf)	% The infinity norm is largest absolute row sum of A.

disp(sprintf('**** Example 2.5 Matrix condition number *'));

A=[.5 1.5 -.5;-.5 2.5 -.5; -.5 -.5 .5];	% Condition number is the ratio of 
			% maximum stretching to maximum shrinking that A produces on a vector.
			% The condition number of a matrix is defined as norm(A)*norm(inv(A)),
norm(A,1)*norm(inv(A),1)	% using norm 1 
norm(A,inf)*norm(inv(A),inf)	% using norm inf 
norm(A,2)*norm(inv(A),2)	% using norm 2 (default in MatLab)
cond(A)				% MatLab condition number of matrix A


disp(sprintf('**** SOLUTION OF A SYSTEM OF LINEAR EQUATIONS *'));

A=rand(3,3)		% 3x3 matrix
b=rand(3,1)		% right hand side
x=A\b			% The MatLab solution      


%%%% Linear least square fitting by quadratic polynomial to points (ti,yi)
t=[-1.0 -0.5 0.0 0.5 1.0];	% five data points
y=[ 1.0  0.5 0.0 0.5 2.0]; 	% measured values 
A=[[ones(1,5)]' t' t'.^2]	% Vandermonde matrix for quadratic polynomial (5x3)
b=y';			% Right-hand side vector is built by measured values.
x=A\b			% The MatLab solution (covers also least squre cases.

disp(sprintf('**** SCRIPTS or FUNCTIONS ****\n'));

%M-files can be either scripts or functions. 
%Scripts are simply files containing a sequence of MATLAB statements. 


%%%%%
%%%Functions make use of their own local variables and accept input arguments.
%A line at the top of a function M-file contains the syntax definition. 
%The name of a function, as defined in the first line of the M-file, should
%be the same as the name of the file without the .m extension. 

%For example, the existence of a file on disk called stat.m with

%    function [mean,stdev] = stat(x)
%    n = length(x);
%    mean = sum(x)/n;
%    stdev = sqrt(sum((x-mean).^2/n));

%defines a new function called stat that calculates the mean and standard 
%deviation of a vector x. The variables within the body of the function are 
%all local variables.
%Subfunctions can be placed also in a single file for a local use.

%%%%%
%%%Control statements:
%if, else, end, break, for, return, switch, while
%%%Syntax 
%    while expression
%        statements 
%    end
%where operation in the expression is: ==, <, >, <=, >=, or ~=.


%%%%%
%%%Nonlinear equations:
x=fzero('x.^3-2*x-5',1)	% Zero of a function of one variable near 1 
xx=fzero('sin(x)',3)	% Zero of sinus function near 3 is Pi
%
p = [1 -6 -72 -27]	% All zeros of polynomial p=1-6*x-72*x^2-27*x^3
r=roots(p)		% Polynomial roots 
%
v=polyval(p,[5 7 -1])	% Polynomial p is evaluated at x= 5,7, and -1 

%%%%%
%%Interpolation:
%%Polynomial
rint=30			% Right side of the observed interval
x = (0: 0.1: rint)';	% Generate a vector of equally spaced x-points (observed interval)
y = sin(x); 		% Evaluate sin(x) at this points (reference curve)
xm = (0: 2: rint)';	% Generate a vector of equally spaced x-points (measurements)
ym = sin(xm); 		% Evaluate sin(x) at this points (measurements)
%
pf = polyfit(xm,ym,6)	% Calculate coefficients of the approximating polynomial of degree 6
f = polyval (pf, x);	% Evaluated values of the approximating polynomial
plot(x,y,'r.',xm,ym,'go',x,f,'k-')	% Plot reference function, measurements 
			% and approximating polynomial.
%interpn		% Multidimensional data interpolation (table lookup) 
%
%%Spline		% Piecewise interpolation
sf = spline(xm,ym);	% Cubic spline interpolation return pp form of spline interpolant
ff = ppval(sf, x);	% Evaluate the approximating spline polynomial
hold on			% Add a new plot in the active window.
plot(x,ff,'b+')
hold off

%%%%%
%% Numerical differentiation:
% Derivative can be calculated by finite difference formulas.
% Alternatively, a polynomial is fitted to function values and explicit derivation 
% of the polynomial is done, to obtain derivatives in other points then knots.
 
%%%%%
%% Numerical Integration:
% Quadrature is a numerical method of finding the area under the graph of a function, 
% that is, computing a definite integral.
il = quad('sin',0,pi)	% Integrate the sine function from 0 to Pi (low order method)

ih = quad8('fun',1,5)	% Integrate the function 'fun' from 1 to 4 (higher order method)
			% The function must be defined in a separate file fun.m as
			% function out1 = fun(x)
			%     out1 = x.^3-2*x-5;   

%%%% Ordinary differential equation solvers (to be used in future Labs)
%ode45, ode23, ode113, ode15s, ode23s, ode23t, 
%ode23tb		% Differential equations (or system of equations) solver                     
%odefile		% Defines a differential equation problem for ODE solvers