CSc
2262 Numerical Methods
Fall 2004
Coates 263, Union 332, CEBA 1302, CEBA Reading
Room, Middleton 109, Middleton 141, Middleton 241, Horticulture 225, New Design
Building 311, and Williams 3rd floor
Three windows – command window (enter
commands and data); graphics window (display plots and graphs), and edit window
(create and edit M-files)
Scalars a = 4; x =
-2.3; A = 0; % case sensitive
variable names
Predefined variables – pi
i ( Ö(-1)
) (can use j for input)
Accuracy format short
(4 decimal places for reals); format long
(15
places)
Arrays 1-dimensional (vector),
2-dimensional (matrix)
a
= [1 2 3 4 5] – row vector
b
= [2;4;6;8;10] or b = [2 4 6 8 10]’ (transpose) – column
vector
c
= [ 1
2
3 ]
Q
= [1 2 3; 4 5 6; 7 8 9] – matrix
Q(2,3)
yields the value ”ans = 6”
who
– lists variables
whos – list variables with
additional information
t = 1:5 yields t = [1 2 3 4 5] %
colon operator
t = 1:0.5:3 yields t = [1.0000
1.5000 2.0000 2.500 3.0000]
t = 10:-1:5 yields t = [10 9 8
7 6 5]
Q(2,:)
yields second row of Q
linspace(r,s,n) yields n points
between r and s
logspace(r,s,n) yields n points
on logarithmic scale between 10r
and
10s
+,
-, *, / are
addition, subtraction, multiplication, and
division
\ is left division for matrices (discussed
later)
-
is negation
^ is exponentiation
Matrix
operators
a
= [1 2 3 4 5]; b = [2;4;6;8;10]; a*b; yields inner
product of 110 while b*a yields a matrix (5 x 5 from outer
product) from matrix multiplication of b (5 x 1) and a (1 x 5). Note that A = [1 2 3; 4 5 6; 7 8 9]; a*A;
yields error message as dimensions are a problem.
A
.^2 yields, thanks to the dot, an element by
element squaring of the elements of A
ones(1:5) yields an array
filled with ones (1)
zeros(0:0.01:5)
yields an array filled with zeros (0)
Built-in
Functions
help log yields help about the
function log
help
elfun yields list of all elementary functions
sqrt, log2, log10,
logm, exp, abs, sin, acos, tanh
adding m to function name at
end yields a matrix operation –
sqrtm
length(t)
yields the number of elements in vector t
size(s)
yields size of matrix s
Graphics
plot (t,v) plots a graph of v
versus t
title(‘Plot of v versus t’);
xlabel(‘Values of t’);ylabel(‘Values of
v’);grid;
colors (blue
b, green g, red r, cyan c, magenta m, yellow y,
black
k)
symbols (point ., circle o,
x-mark x, plus +, star *, square s,
diamond
d, triangle
down
v, triangle up ^, triangle left <, triangle right >,
pentogram p,
hexagram
h)
line types (solid -, dotted :,
dashdot -., dashed --)
Help
lookfor algorithm yields all
references in online help with word
“algorithm”
Files
M-file – file containing MATLAP
“program”
Edit, enter program, save
Debug
Run
– name file in command
Script file
Function file
function
outvar = functionname(argumentlist)
statements
outvar
=value
save in file named functionname
I/O
x = (‘promptstring’)
disp(value)~
fprintf(format like C)
x
= [1 2 3 4 5]; y = [20.4, 12.6 17.8 88.7 120.4]; z = [x;y];
fprintf(‘ x
y\n’); fprintf(‘%5d, %10.3f\n’ z);
Structured
Programming
if
conditional statement(s) end
if
conditional statement(s) else
statement(s) end
if
conditional statement(s) elseif
conditional statement(s) elseif
conditional
statement(s) … elseif statement(s) else
statement(s)
end
conditionals
== (equals), ~= (not equals), <
(less than), > (greater than),
<=
(less than or equal), >= (greater
than or equal)
Loops
for index – start; step; finish
statement(s) end
default step is 1, e.g. for i = 1:5 disp(i) end
while condition statement(s) end
while (1) statement(s) if condition break, end statement(s)
end
Nesting
function quad = quadroots(a,b,c);
if a == 0
if b ~= 0 r1 = -c/b else error(‘trivial solution, try
again’)
end
else
d = b^2 –
4*a*c;
if d >= 0 r1 = (-b + sqrt(d))/(2*a); r2 = (-b –
sqrt(d))/(2*a);
else r1 = -b/(2*a);
r2 = r1; i1 = sqrt(abs(d))/2*a; i2 = -
i1
end
Passing
Functions as Parameters
Outvar = feval(‘funcname’, arglist)
function
f = func(x)
f
= 0.125*x.^3 – 1.125*x.^2 + 2.75*x + 1;
function
favg = funcavg(a,b,n)
x
= linspace(a,b,n); y = feval(‘func’), x); favg = mean(y);
Example
Consider
someone who wants to do bungee jumping. To calculate the velocity over time,
considering both gravity and drag (resistance due to wind), we have dv/dt = g –
cd v2/m,
where v = velocity (m/s), t = time (s), g = acceleration due to gravity = 9.81
m/s2), cd
= second-order drag coefficient (kg/m), and m = jumper’s mass (kg). This can be
solved to yield v(t) = Ö(gm/cd)
tanh(Ö(gm/cd)t),
where tanh is the hyperbolic tangent function, i.e., tanh(x) = (ex–e-x)/
(ex+e-x). But, we can approximate this analytical or closed-form solution by numerical methods using vi+1 = vi_+dvi/dt
Dt, with Dt = ti+1-ti, i.e., Euler’s method (we will discuss this
method in detail later). Note that dvi/dt can be approximated by (vi+1 -vi)/(ti+1-ti)
= g – cdv(ti2)/m.
function dv =
derive(t,v,m,cd)
g =
9.81;
dv
= g–(cd/m)*v^2;
function vend = velocity2(dt, ti, tf, vi, m, cd)
t =
t1; v = vi; h = dt;
while
(1)
if t+dt > tf, h = tf –t; end
dvdt
= derive(t,v,m,cd);
v
= v+ dvdt*h;
t
= t+h;
if t >= tf, break, end
end
vend
= v
One might try this
with m = 68.1, cd = 0.25, ti = 0, tf = 10, and
dt
= (tf-ti)/n with n=5, with vi = 0.