### old simple class notes

```math403
Introductory MATLAB
Zheng Chen
Spring 2010
1
Environment of matlab
•
•
•
•
3 windows:
Command window
History window
Editor window
2
Executable file form
• Files of form .m for functions and main
programs
• All files to solve one problem should be in one
directory
• Pwd, cd, ls, cd.., cp, mkdir, rmdir,
• Help while
3
Matlab as a calculator
•
•
•
•
•
•
2+3/4.0*5+(-2)^3+ log(3.4)
Sin(pi/4);
1.2+3i
Format of numbers
Format long
Format short
4
variables
•
•
•
•
•
X=3-2^3;
Assign the value to x
Y=x+4; assign the value to y
Pi=3.1415…, eps=2.2204e-16
I, j have the value sqrt(-1);
5
Built in functions
•
•
•
•
Trig
Log
Sqrt(x);
exp
6
vectors
•
•
•
•
•
•
•
•
Row vector: V=[1,2 3]
U=[-5, 6, 12]
V+/-U;
W=[7,8];
Cd=[2*w, 3V];
V(2); it is the same as V(1,2);
Column vector w=[7;8;9]
W(2) =w(2,1);
7
Colon notations
• V=1:8;
• V=3:7;
• V=.32:0.1:-2
• Extracting
• A=[1,2 ,3 ,4,5 ,6 ,7 ,8]
• What is the b=A(4:7) and c=A(1:3:7)
8
transpose
• Row vector to column vector or vice verse
• General matrices with size m*n
• W’ is of size n*m if W has size m*n
9
Scalar product of vectors u*v
•
•
•
•
•
Scalar product:
u=[a, b, c], v=[l, m,n]’
In matlab, u*v=al+bm+cn; it is a number
sqrt(u*u’) ; and norm(u);
Recall what is the meaning of norm?
10
Dot product .*
• u=[a, b, c], v=[l, m, n]
• u.*v is a vector with the same size as u and v;
• simply multiply the corresponding elements of
two vectors
• u.*v named as dot product
11
Dot Division ./
• ./ is defined to give element by element
division
• u=[a, b, c], v=[l, m, n];
• u./v=[a/l , b/m, c/n];
• Try 1./[2, 3, 4];
12
Dot Power of Arrays (.^)
• u = [ 10, -11, 12];
• u.^2;
• .* ; ./ ; .^ mean component wise
operations
13
Functions applied to arrays
• u = [ 10, 11, 12];
• sin(u), sqrt(u);
14
Matrices –
Two-Dimensional Arrays
• Input: a=[1 2 3; 5 6 7];
• or a= [1 2 3
5 6 7];
b=[1 2; 3 4; 5 6];
• size(a) = dimensions of matrix, that is, number
of rows and number of columns
• size(a) is itself one vector with size 1, 2.
15
Special Matrices
• ones(2,3); zeros(2,3); eye(3);
• d = [-3 4 2], D = diag(d);
• given F = [0 1 8 7; 3 -2 -4 2; 4 2 1 1],
find diag(F);
• c=[0 1; 3 -2; 4 2]; x=[8;-4;1];
find g = [c x] (extra column)
find h=[c; k] (stacked c and k) where k=[7 8]
• Try spy(F), grid for fun
16
Extracting Bits of Matrices
• J =[ 1 2 3 4
5 6 7 8
9 0 3 11];
Find J(2,4); J(:,2); J(3,:); J(1:2, 2:4);
17
Dot product of matrices (.*)
• The dot product works as for vectors:
corresponding elements are multiplied
together
18
Matrix{vector products
• An (m x n) matrix times an (n x k) matrix is a
(m x k) matrix.
• A = [5 7 9; 1 -3 -7]; x = [8; -4; 1]; find A*x
• Practice: choose two matrices A and B, find
A*B, (A*B)’ and B’*A’
• Practice: choose matrices A, B and C, find
A*(B*C) and (A*B)*C
19
Sparse matrices
•
•
•
•
•
•
•
•
>> i = [1, 3, 5]; j = [2,3,4];
>> v = [10 11 12];
>> S = sparse (i,j,v); T=full(S);
save memory, only save diagonal elements
>> n = 5;
>> l = -(2:n+1)'; d = (1:n )'; u = ((n+1):-1:2)';
>> B = spdiags([l d u],-1:1,n,n);
>> full(B) (%l,d,u are column vectors)
20
System of linear equations
• Ax=b
• x = A \ b; (“\backslash")
• try inv(A) and inv(A)*b
21
Overdetermined system
of linear equations
• generally speaking, Ax=b has no solution
when # of equations > # of unknowns
• r = Ax-b;
• Least square solution: make r minimum
• normal equations:
Cx = d; where C = A’A and d = A’b
• matlab gives the solution in its routine library
x = A\b
22
• A =[0.0000 0.0000; 0.0001 0.0000;
0.0013 0.0000; 0.0030 0.0000; 0.0075 0.0001]
• b = [5 50 500 1000 2000];
• X=A\b’
23
For Loop
• for i=1:10
disp([i i^2 i^3])
end
• for i=1:4
for j=1:4
disp([i j i^2+j^2])
end
end
24
example
• for i=1:4
for j=i:4
disp([i j i^2+j^2])
end
end
25
example
• sum=0.0; bigsum=0.0;
• for i=i:1000
bigsum =bigsum+1/i;
end
for i=i:1000
sum =sum+1/i^2;
end
disp([ bigsum sum])
26
Questions
• What is the sum of 1+1/2+1/3+… until the
term where 1/n<=10^(-3) ?
• How many terms we have in the expansion of
1+1/2+1/3+… if the partial sum is just >=2 ?
27
While loop
• % format long
i = 1; % start the counter at 1
total = 0; % initialize the sum
while(1/i >= 1e-3)
% condition
total = total + 1/i;
% update the sum
i=i+1;
end
disp([i 1/i total])
% display the number of
loops, 1/i, and the current sum
28
While loop(condition changed)
i = 1; % start the counter at 1
total = 0; % initialize the sum
while(total <= 2) % condition changed
total = total + 1/i
% update the sum
i= i+1
end
% disp([I total])
disp([i-2 total-1/(i-1)])
% display the
number of terms, and the final sum
29
Example to use while loop
• Problem: what is the greatest value of n that can
be used in the sum
1+(1/2)^2+ (½)^2 +… + (1/n)^2
and get a value of less than 2?
S = 1; n = 1;
while S+ (1/(n+1))^2 < 2
n = n+1; S = S + (1/n)^2;
end
• [n, S]
30
Logicals
• true = 1, false = 0
• If x<1.0
elseif
31
if...then...else...end
• if logical test 1
commands to be executed if test 1 is true
elseif logical test 2
commands to be executed if test 2 is true
but test 1 is false
...
end
32
• R = input('What is your name:','s')
if R=='cuana'
else
'bye'
end
33
• Problem: find out all triples of integers, at
most 20, as the 3 sides of one right triangle,
like 3, 4, 5;
34
Using for and if loop
for i=1:20
for j=i:20
for k=j:20
if k^2==i^2+j^2
[i, j, k]
end
end
end
end
35
Functions defined by yourself
•
•
•
•
•
% Compute the area of a triangle whose
% sides have length a, b and c.
% Inputs: a,b,c: Lengths of sides
% Output: A: area of triangle
% Usage: Area = area(2,3,4);
function [A] = area(a,b,c)
s = (a+b+c)/2;
A = sqrt(s*(s-a)*(s-b)*(s-c));
36
• Problem: find the are of one triangle with
given 3 numbers, you need to decide whether
these three numbers can be the 3 sides of one
triangle.
for example, 3, 4 and 5 are 3 sides of one
triangle, 3, 4 and 8 are not 3 sides of one
triangle.
given a<=b<=c , criteria for them to be 3 sides
of a triangle: c<a+b & a> c-b.
37
Functions countinue
• function [A] =triValue(a,b,c)
s = (a+b+c)/2;
A = sqrt(s*(s-a)*(s-b)*(s-c));
• save the file as triValue.m
• Another function: two results are provided
function [cir, A] =triValue(a,b,c)
s = (a+b+c)/2;
cir=s*2;
A = sqrt(s*(s-a)*(s-b)*(s-c));
38
• function area=tri_area(a, b, c)
dd=min([a,b,c]);
aa=max([a,b,c]);
b=a+b+c-dd-aa;
a=aa; c=dd;
if c>(a-b) & b+c>a
'the given numbers are 3 sides of one triangle, and the area is: '
s=(a+b+c)/2;
area=sqrt(s*(s-a)*(s-b)*(s-c))
else
'the given numbers are not 3 sides of one triangle'
end
39
Another example of functions
• Fibonnaci sequence is dened by
f1  0; f 2  0; f n  f n 1  f n  2
n=3,4, …
input: integrer
output: f n
40
Method 1: File Fib1.m
function f = Fib1(n)
% Returns the nth number in the
% Fibonacci sequence.
F=zeros(1,n+1);
F(2) = 1;
for i = 3:n+1
F(i) = F(i-1) + F(i-2);
end
f = F(n);
41
function f = Fib2(n)
% Returns the nth number in the
Fibonacci sequence.
if n==1
f = 0;
elseif n==2
f = 1;
else
f1 = 0; f2 = 1;
for i = 2:n-1
f = f1 + f2;
f1=f2; f2 = f;
end
end
42
Recursive programming
function f = Fib3(n)
% Returns the nth number in the
% Fibonacci sequence.
if n==1
f = 0;
elseif n==2
f = 1;
else
f = Fib3(n-1) + Fib3(n-2);
end
43
% The final version uses matrix powers. The vector y has
two components,
 fn 
y  

 f n 1 
% Returns the nth number in the Fibonacci sequence.
function f = Fib4(n)
A = [0 1;1 1];
y = A^n*[1;0];
f=y(1);
44
Which method you like most?
• Which algorithm is most efficient?
• method 3(recursive programming) no efficient
45
More Built-in Functions
•
•
•
•
•
x = pi*(-1:3), round(x)
Try fix(x); floor(x); ceil(x); sign(x), rem(x,3)
A = [1:3; 4:6; 7:9];
s = sum(A), ss = sum(sum(A));
x = [1.3 -2.4 0 2.3], max(x), max(max(A));
min(x);
• y = rand, Y = rand(2,3)
46
Use rand
% simulates "n" throws of a pair of dice
% Input: n, the number of throws
% Output: an n times 2 matrix, each row
% referring to one throw.
% Useage: T = dice(3)
function [d] = dice(n)
d = floor(1 + 6*rand(n,2));
%% end of dice
47
Function ‘find’
• A = [ -2 3 4 4; 0 5 -1 6; 6 8 0 1];
• k = find(A==0);
To interpret the result we have to recognize
That “find" first reshapes A into a column vector;
this is equivalent to numbering the elements of A
by columns as in
1 4 7 10
2 5 8 11
3 6 9 12
Now try n = find(A <= 0); and A(n)
48
• A=[1 4 7 10
2 5 8 11
3 6 9 12];
n = find(A <= 0);
Try A(n);
Try m = find( A' == 0);
49
plot
•
•
•
•
•
linspace (a,b,n) producing n+1 points;
N = 10; h = 1/N; x = 0:h:1; doing the same.
X=linspace(0,1,110);
y = sin(3*pi*x);
Plot(x,y)
50
Plot with title
• Using string inside of ‘ ’
title('Graph of y = sin(3pi x)')
xlabel('x axis')
ylabel('y-axis')
grid
off (% off grid)
51
Plot-line styles & color
• plot(x,y,'w-');
• Colors
y
yellow .
m
magenta
c
cyan
R
red
g
green
b
blue
w
white
k
black
line styles
point
o
circle
x
x-mark
+
plus
solid
*
star
:
dotted
-.
dashdot
-dashed
• Try plot(x,y,'b*‘)
52
• plot(x,y,'w-'), hold on
• plot(x,y,'gx'), hold off
• save or print the graphs in .eps form
• clf, N = 100; h = 1/N; x = 0:h:1;
• y = sin(3*pi*x); plot(x,y)
• axis([-0.5 1.5 -1.2 1.2]), grid
53
subplot
for n = 1:8
subplot(2,4,n), plot(x,sin(n*pi*x))
end
54
Characters, Strings and Text
• %The ability to process text in numerical
processing is useful for the input and output
of data to the screen or to disk-les.
• The assignment, t1 = ‘abc ‘, t2 = ‘efghijk‘;
• t3 = [t1,t2];
• T3 is abcdefghijk;
55
Conversion among strings,
integers and real numbers
• 'str2num‘, 'int2tr' and 'num2str‘;
• Example:
N = 11:20; ncube = N.^2;
['The perfect squares of integere from ',
int2str(ncube(1)),' to ', int2str(ncube(10)), '
are in the following: ']
num2str(ncube)
% fprintf(num2str(ncube))
% disp(ncube)
56
Plotting Surfaces
• Plot the surface of f(x; y) = (x -3)^2-(y-2)^2
for 2<= x <= 4 and 1<= y<= 3.
[X,Y] = meshgrid(2:.5:4, 1:.5:3);
Z = (X-3).^2-(Y-2).^2;
mesh(X,Y,Z)
57
f=-xyexp(-2(x^2+y^2)), -2<=x,y<=2;
[X,Y] = meshgrid(-2:.1:2,-2:.2:2);
f = -X.*Y.*exp(-2*(X.^2+Y.^2));
figure (1)
mesh(X,Y,f), xlabel('x'), ylabel('y'), grid
figure (2), contour(X,Y,f)
xlabel('x'), ylabel('y'), grid, hold on
58
To locate the maxima of the
“f” values on the grid
fmax = max(max(f))
kmax = find(f==fmax)
Pos = [X(kmax), Y(kmax)]
plot(X(kmax),Y(kmax),'*')
text(X(kmax),Y(kmax),' Maximum‘)
% the ‘find’ above gives the index numbers of
matrices X and Y
59
timing
% Matlab allows the timing of sections of code
by providing the functions tic and toc.
tic
x=0;
for j=1:1000
x=x+j^3;
end
toc
60
Data Files
In a file 'table.dat' , data is as following:
1 56
2 58
3 65
4 67
fid = fopen('table.dat','r');
% action is reading, file identfier fid
a = fscanf(fid,'%1d%2d');
% one with 1 digits and one with 2 digits
fclose(fid); % close the file with file identier fid.
A = reshape(a,2,5)';. % change the size of a to 2times5
61
HW and projects
• Pob1: using the mid point method to find the
approximate solution of any equation in one
interval for example: 4x+3sin(pi/2*x)-5=0 in
[0,1]
• Pob2: using Newton’s method to find the
approximate solution of any equation in one
interval for example: 4x+3sin(pi/2*x)-5=0 in
[0,1]
62
• Pob3: find the determinant value of any given
people in input the matrix row by row)
• Pob4: using echelon method to find the solution
of any linear system with 2 unknowns; your
program should ask people to input the
coefficients of the system and the program
should be valid to any such system, in other
words, your program should contain 3 possible
situations of solutions.
63
• Given 3 points, find the distance from the
third point to the line though the other two