Transcript M ATLAB Matrix Laboratory

MATLAB
Matrix Laboratory
Internet Resources:
Dr. H-C Chen’s CVEN 302 Web Page
http://ceprofs.tamu.edu/hchen/cven302.html
For MATLAB, select chap01b.ppt, or click
http://ceprofs.tamu.edu/hchen/cven302/chap01b.ppt
Programming
The ideal style of programming is structured
• Break down a larger goal into tasks
• Develop a module for each task
• A module has a single entrance and exit
• Modules are reusable
Use MATLAB to solve a problem
•
•
•
•
•
•
State the problem clearly
Describe the Input/Output (I/O)
Algorithm - Numerical Method
Develop a MATLAB program
Debugging and Testing
Documentation
If needed, work the problem by hand
Compiler and Program Execution
input
Computer
Program
Compiler
Machine
language
program
Link/
Load
Execute
output

  
compilation
linking/loading
execution
MATLAB Execution
Interactive environment - does not require formal
compilation, linking/loading, and execution
Scripts - develop and execute m-files that contain
MATLAB commands
Important!
You must include comments in your
programs you turn in - otherwise we will
have great difficulty knowing what you
are thinking and doing. You will have
difficulty knowing what you were doing
if you need to reuse the programs.
MATLAB’s 3 Windows
Command Window
- enter commands and data
- print results
Graphics Window
- display plots and graphs
Edit Window
- create and modify m-files
Two simple MATLAB programs (also
called scripts)
Filename: lec2a.m
% Find sum of 1 to b
% Comments begin with “%”
% cd y:\cven302\
clear all
sum = 0;
b = input('what is the value of b?
for i = 1 : b
sum = sum + i
end
Need extension .m
called m-file
')
Filename: lec2b.m
% plot 1/x v.s. tanh x
clc
a = linspace(0,5);
f = 1./a;
g = tanh(a);
plot(a,f,'-', a, g,'--')
axis([0 3 0 2])
xlabel ('x')
ylabel ('y')
title ('Plot 1/x and tanh x')
text(0.7, 1.6, '1/x')
text(2, 1.1, 'tanh x')
Now we have seen Command Window,
Graphics Window, and Edit Window
There is another window which is very useful:
the Help Window
You can also type help in Commend Window
e.g.,
help input
» help input
INPUT Prompt for user input.
R = INPUT('How many apples') gives the user the prompt in the
text string and then waits for input from the keyboard.
The input can be any MATLAB expression, which is evaluated,
using the variables in the current workspace, and the result
returned in R. If the user presses the return key without
entering anything, INPUT returns an empty matrix.
R = INPUT('What is your name','s') gives the prompt in the text
string and waits for character string input. The typed input
is not evaluated; the characters are simply returned as a
MATLAB string.
See also KEYBOARD.
MATLAB’s basic component is a matrix
i.e., all variables are treated as matrices
All operations are vectorized (optimized
for vector use)
Loops run much slower in MATLAB
than in Fortran (not a vector operation)
Scalar, vector, and matrix
Matrix:
1 2 3 4 
a   3 5 7 9


2 4 6 8
m  n matrix
m rows (m = 3)
n columns (n = 4)
Scalar: k = 5
11 matrix
Vector:
b  1 2 3 4
14 matrix, or
a row vector
1 
c   3
 
2
31 matrix, or a
column vector
Variable
• may consist up to 31 characters
• starting with a letter and followed by any
combination of letters, digits, and underscores
(Variable names must start with a letter)
• punctuation marks and spaces should not be
included
e.g., CVEN_302, TIME, Time, time, velocity, force_x,
force_y
Data Type
• All numbers are double precision
• Text is stored as arrays of characters
• You don’t have to declare the type of data
(defined when running)
• Variables are case sensitive
e.g., CVEN_302, TIME, Time, time, velocity, force_x,
force_y
Arithmetic operation of scalars
Addition:
a+b
multiplication: a*b
exponentiation: a^b= ab
subtraction: a-b
division: a/b
Elementary math functions
abs(x) : absolute value
sqrt(x) : square root = x^(1/2)
sin(x) : sine
asin(x) : inverse sine
sinh(x) : hyperbolic sine
asinh(x) : inverse hyperbolic sine
log(x) : natural logarithm
log10(x) : common logarithm
exp(x) : exponental: ex (e = 2.7183….)
» a = 5;
» b = 3;
» a*b
ans =
15
» d = ans * 3
d=
45
» e = ans * 3;
»e
e=
45
» a^2
ans =
25
» a/b
ans =
1.6667
» format long
» a/b
ans =
1.66666666666667
» format short
» a/b
ans =
1.6667
» exp(2)
ans =
e2
7.3891
» 10^2
ans =
100
» log(100)
ln 100
ans =
4.6052
» log10(100)
ans =
log 100
2
» pi
ans =
3.1416
; - result not displayed
format long - 15 digits
format short - 5 digits
ans - result from previous calculation
pi - 
Ctrl-c : terminate a running program
Array Operations
• An array operation is performed element-byelement - Need “.” in front of the operator
MATLAB: C = A.B;
C(1)  A(1)  B(1);
C(2)  A(2)  B(2);
C(3)  A(3)  B(3);
C(4)  A(4)  B(4);
C(5)  A(5)  B(5);
Arithmetic operation of Arrays
Addition:
a+b
multiplication: a.*b
exponentiation: a.^b
subtraction: a-b
division: a./b
= ab
MATLAB variables and matrix operation
b  4 5 6
» b = [4 5 6]
b=
4 5 6
» c = [1
3
2]
c=
1
3
2
1 
c   3
 
2
» c = [1; 3; 2]
c=
1
; - a new line
3
2
» d = c'
 - transpose
d=
1 3 2
b  4 5 6
1 
c   3
 
2
»b*c
ans =
d  1 3 2
31
» b.*d
ans =
b.*d = [4*1 5*3 6*2]
4 15 12
b./d = [4/1 5/3 6/2]
» b./d
ans =
4.0000 1.6667 3.0000
» f(2,3)
» f = [b; d; [3 6 8]]
ans =
f=
2
b
4 5 6
» size(f)
d
1 3 2
[3 6 8]
3 6 8
size - size of a matrix
Matrix Operation (conti.)
x  1 2 3 ;
y  7 9 4
z  x y  1 2 3 7 9 4
u  x ; y   1 2 3
7 9 4
s = [1 2 3 4; 5 6 7 8]
t = [1 2 3 4....
5 6 7 8]
continue
Colon Operator
Creating new matrices from an existing matrix
C = [1,2,5; -1,0,1; 3,2,-1; 0,1,4]
F = C(:, 2:3) = [2,5; 0,1; 2,-1; 1,4]
all columns
E = C(2:3,:) = [-1 0 1; 3 2 -1]
G = C(3:4,1:2) = [3,2; 0,1]
 1
 1
C
 3

 0
5
0
1

2  1

1 4
2
5
2
0
1

F
2  1


1
4


rows 2 to 3
1
 1 0
E

3
2

1


 3 2
G

0
1


a : step : b
gives number from a to b with
the specified step between elements
step = 1 if not specified
» x = 1:10
x=
1 2 3 4
» t = 1 : 2 : 10
t=
1 3 5 7
» k = 5:-1:-3
k=
5 4 3 2
5
6
7
8
9
0
-1
-2 -3
9
1
10
What if the step is not easy to calculate, or is an odd number?
y = linspace (a, b, n_pts)
» y = linspace (0, 1, 4)
y=
0 0.3333 0.6667
2 end points are
included in n_pts.
1.0000
y = linspace (0,1)
a number of 100 for
n_pts is assigned.
Some useful special matrices
» eye (3)
ans =
1 0 0
0 1 0
0 0 1
» ones (3)
ans =
1 1 1
1 1 1
1 1 1
» ones (1,4)
ans =
1 1 1
eye (m,n)
ones (m,n)
zeros (m,n)
1
» zeros (2,3)
ans =
0 0 0
0 0 0
eye (n)
ones (n)
zeros (n)
Plotting
plot (x, y)
plot(x1, y1, x2, y2)
plot (x, y, ‘color symbol line style’)
» x = linspace(0, 2*pi);
» y = sin (x);
» z = cos (x);
» plot (x, y)
figure or figure (#) : open a figure
» plot (x, y, x, z)
» figure (2)
» plot (x, y, 'r o -'); grid on add grids
try help plot
» hold on
red, circle, solid
» plot (x, z, 'b x :')
blue, x-mark, dotted
try grid off
Plotting (continue)
xlabel (‘ label ‘)
ylabel (‘ label ‘)
title (‘ title of the plot ‘)
text ( x_location, y_location, ‘ text ‘)
axis ([ x_min x_max y_min y_max ])
» xlabel ('x')
» ylabel ('y')
» title ('sine and cosine')
» text (2, 1, 'This is a sine fuction')
» axis ([0 2*pi -2 2])
» hold off
  - text string
hold on: hold the
plots to avoid
overwriting
Plotting (continue)
subplot ( m, n, figure number ) - break the
figure window into m by n small figures, and
plot the specified figure
» figure
» subplot (3, 2, 1)
» plot (x, y)
» subplot (3, 2, 2)
» plot (x, z)
» subplot (3, 2, 4)
» plot (x, y-z)
1
3
5
2
4
6
semilogx (x, y) logarithmic scales for the x axis.
semilogy (x, y) logarithmic scales for the y axis.
loglog (x, y) logarithmic scales for the x and y axes
x = 10.^[-1:0.1:2];
y = exp(x);
figure(1)
semilogy(x,y,':^')
figure(2)
loglog(x,y,'-s')
grid on
More control on plotting
get(H) get object properties
set(H,'PropertyName',PropertyValue,...)
set object properties
figure(2)
h = loglog(x,y,'-sr')
grid on
get(h)
set(h, 'LineWidth',2)
set(h, 'MarkerFaceColor', [0 0 1])
Some useful commands
who
list of current variables
clc
clear the command window
clf
clear the graphical windows
clear x
clear the variable x
clear all clear all variables
close
close the current figure
close all close all figures
cd y:\cven302\
change directory
dir
list all files
what
list all m-files
CTRL-C Abort
semicolons (;) at the end of line: Calculation results
will not be displaced on the command window
Initializing Variables
•
•
•
•
Explicitly list the values 
reads from a data file
uses the colon “:” operator 
reads from the keyboard
Input and output (I/O)
a = input (‘ enter a value or string ‘)
- wait for input from the keyboard
load file_name
- input an existing data file (ASCII format)
diary file_name
- save anything in the Command Window
diary off
save file_name variable(s) -ascii
- save the data of the variable(s) to a file in
ASCII format
You may need to go to certain directory before loading and/or
saving files
» diary tmp.dat
» a = [ 1 2 3; 4 5 6]
a=
1 2 3
4 5 6
» save a_file.dat a -ascii
» load a_file.dat
» b = a_file
b=
1 2 3
4 5 6
» d = input ('give a value to d: ')
give a value to d: 5
d=
5
» diary off
Save in the
tmp.dat file
M-Files: Scripts and Functions
• You can create and save code in text files using
MATLAB Editor/Debugger or other text editors
(called m-files since the ending must be .m)
• M-file is an ASCII text file similar to FORTRAN or
C source codes ( computer programs)
• A script can be executed by typing the file name, or
using the “run” command
Functions
distinguished from the script that the first line is of the
form
function x = function_name (input arguments)
function [x, y] = function_name (input arguments)
A function has to be stored as a stand-along file ended
with “.m”. The name of the function is usually (but
not necessary) the same as the name of the file.
Difference between scripts and functions
Scripts share variables with the main workspace
Functions do not
Approximate integral f(x) = 1/x3 using basic trapezoid rule
File name: trap_ex.m
function t = trap_ex(a, b)
t = (b - a) * (a^(-3) + b^(-3)) / 2;
» y = trap_ex (1, 3)
y=
1.0370
»
y
b
a
1
dx
3
x

b  a
 f a   f b 
I
2
Trapezoidal rule
f(b)
f(a)
I  width  average height
a
b
Approximate integral f(x) = 1/x3 using basic trapezoid rule
with n steps used
File name: trap_ex2.m
function t = trap_ex2(a, b, step)
dx = (b - a)/step;
x1 = a;
t = 0;
for i = 1 : step
x1 = a + (i-1)*dx;
x2 = a + i*dx;
t = t + (x2 - x1) * (x1^(-3) + x2^(-3)) / 2;
end
» y = trap_ex2 (1,3,1)
y=
1.0370
» y = trap_ex2 (1,3,2)
y=
0.6435
» y = trap_ex2 (1,3,100)
y=
0.4445
» y = trap_ex2 (1,3,1000)
y=
» x1
Why?
0.4444
??? Undefined function or variable 'x1'.
» format long
» y = trap_ex2 (1,3,1000)
y=
Exact solution = 0.444444444444 ….
0.44444543209743
feval - evaluate function specified by string
my_func.m
function y = my_func(x)
% function 1/x^3
y = x.^(-3);
basic_trap.m
function q = basic_trap(f, a, b)
% trapezoid rule with input function f
ya = feval(f, a);
yb = feval(f, b);
q = (b - a)* (ya + yb)/2;
» feval('my_func', 1)
ans =
1
» y = basic_trap ('my_func', 1,3)
y=
1.0370
Decision Making: Control Flow
(a)
For Loops
(b)
While Loops
(c)
If-Else Structures
For Loops
for x = array
commands
end
z = 0;
for i = 1:10
y(11-i) = i;
z = z+i;
end
»y
y=
10 9
» disp(z)
55
8
7
Without this line you will get:
Warning: Reference to uninitialized
variable z.
6
5
4
3
2
1
disp (x) - display the value of x
Filename: lec2a.m
% Find sum of 1 to b
% Comments begin with “%”
% cd y:\cven302\
clear all
sum = 0;
b = input('what is the value of b?
for i = 1 : b
sum = sum + i
end
')
While Loops
while expression (is true)
commands
end
eps = 1;
count = 0;
while (1+eps) > 1
eps = eps/2;
count = count + 1;
end
eps = eps*2
display(count-1)
eps =
2.2204e-016
ans =
52
Determining machine
epsilon
Floating point relative accuracy
52 digits in double precision (64 bit)
Floating point
exponent
N  .d1d2d3d4    d p Be
a number
digit
base
d1 0: normalized floating-point system
Single precision: 32 bits (23 for digit, 8 for exponent, 1
for sign
Double precision: 64 bits (52, 11, 1)
eps = 1;
count = 0;
while (1+eps) > 1
eps = eps/2;
count = count + 1;
end
eps = eps*2
display(count-1)
What if we omit this operation?
The loop becomes an “infinite
loop” because the statement is
always true.
eps = 1;
count = 0;
while (1+eps) > 1 & count < 100
eps = eps/2;
count = count + 1;
end
eps = eps*2
display(count)
Use “break” to terminate
while and for loops
prematurely (or set a limit
in while loops)
If-Else Structures
if expression (is true)
commands
end
if expression
commands 1
else
commands 2
end
if expression 1
commands 1
elseif expression 2
commands 2
elseif expression 3
commands 3
:
:
else
commands n
end
Relational Operators
<
<=
>
>=
==
~=
less than
less than or equal to
greater than
greater than or equal to
equal to
not equal to
Logical Operators
&
|
~
and
or
not
eps = 1;
count = 0;
while (1+eps) > 1
eps = eps/2;
count = count + 1;
if count > 100
break
end
end
eps = eps*2
display(count)
eps = 1;
count = 0;
while (1+eps) > 1 & count < 100
eps = eps/2;
count = count + 1;
end
eps = eps*2
display(count)
a = input('what is the value of input data?
if a > 0
sign = 1;
elseif a < 0
sign = -1;
else
sign = 0;
end
disp(sign)
')
Print figures to image files
print - print out current figure
print -djpeg filename - save current figure in
jpeg format
print -djpeg# filename - save current figure in
jpeg format with # (resolution level)
between 0 and 100 (default 75)
print -dtiff filename - save current figure in
tiff format
print -dpsc filename - save current figure in
color PostScript format
See help print for more
Filename: lec2b.m
% plot 1/x v.s. tanh x
clc
a = linspace(0,5);
f = 1./a;
g = tanh(a);
plot(a,f,'-', a, g,'--')
axis([0 3 0 2])
xlabel ('x')
ylabel ('y')
title ('Plot 1/x and tanh x')
text(0.7, 1.6, '1/x')
text(2, 1.1, 'tanh x')
» cd y:\cven302
» print -djpeg myfigure1
» print -djpeg100 myfigure2