CSc 2262 Numerical Methods

Fall 2004

 

MATLAB Tutorial

 
Computer Laboratories on Campus with Matlab

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

Environment

          Three windows – command window (enter commands and data); graphics window (display plots and graphs), and edit window (create and edit M-files)

Assignment

          Scalars                 a = 4; x = -2.3; A = 0;          % case sensitive variable names

          Predefined variablespi

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”

More

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

Mathematical operators

                    +, -, *, /  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.