MATH2071: LAB #2: MATLAB

TABLE OF CONTENTS

• Introduction
• Messing with MATLAB
• MATLAB Data
• Variables
• Variables Are Vectors
• Vector Operations
• Polynomials
• Script Files
• Function Files
• Graphics
• ASSIGNMENT

Introduction

If you already know this stuff, skip ahead to Simple ODE Methods.

In the last lab, we started using MATLAB without even thinking about what we were doing. It wasn't too hard to get started, but now we need to learn a few facts and rules, so that we can

• understand the data structures
• generate and manipulate vectors,
• write reusable function routines,
• make some simple graphics.

Messing with MATLAB

MATLAB is very good at evaluating a formula, or setting up a little algorithm. Let's see what it's like to use MATLAB. Start the program by typing matlab (...or octave). Determine the hypotenuse of the (3,4,5) triangle:

x = 3;
y = 4;
z = sqrt ( x^2 + y^2 )

x = 1;
f = 1 / ( 1 + x^2 )
x = 4;
f

Now here's an interesting question. If we start out with the number x1 = 1, and then invert it and add it to 1 to get x2, and so on, what is the limit? We could try this:

x1 = 1
x2 = 1 + 1 / x1
x3 = 1 + 1 / x2
x4 = 1 + 1 / x3

x = 1
for i = 1 : 10
x = 1 + 1 / x
end

If you want to exit MATLAB now, type quit, although you can just leave it running while we continue with more exercises. But now let's pause for a minute and discuss some details.

MATLAB Data

We never thought about what goes on when we say

x = 1

Exercise: look how the same variable can hold different types of data:

x = 1
x = 1 / 2
x = exp ( 1 )
x = sqrt ( - 1 )
x = 2 + sqrt ( - 1 )
x = 1 / 0
x = 0 / 0
x = 'Howdy!'
x = [ 1, 2, 3 ]
x = [ 1, 2, 4; 4, 5, 6 ]

What MATLAB actually does is to store everything as a double precision complex array. This means that the value stored in x has about 14 digits of accuracy, that it can be set to a complex value if we like, and that x can be adjusted to accommodate a string or table of values instead of just one value.

Exercise: Verify that a MATLAB variable has about 14 digits of accuracy by computing the machine epsilon, also called the machine accuracy or unit round. The machine epsilon is the smallest positive number that can be added to 1 and produce a result larger than 1. Assume that epsilon is a power of 2, and initialize it to 1. What do you have to test to see whether to divide by 2?

         epsilon = 1;
         while ( 1 + epsilon > 1 )
           epsilon = epsilon / 2;
         end
         epsilon
      

You can test your value of epsilon in two ways:

• type the command eps, which will print out the built-in value of the machine epsilon;
• Compute the value of the formula

  ( 1 + epsilon ) - 1

Variables

One of the simplest commands in MATLAB has the form:

x

y
n
banana
f(x)
max ( x, y )

To see the current value of a named variable, just type its name.

Of course, this raises the question of how to put values into variables. MATLAB uses the = sign to do this, which is a standard computer science usage. The variable on the left hand side of the equals sign is assigned the value of the constant, variable, or expression on the right hand side. It is common for a variable name to appear on both sides of the assignment statement.

x = 1
y = x + 2
x = x + 1
z = cos ( x )

You can name your variables anything you want to, but you should be careful not to use certain names that MATLAB has already defined. These include:

• ans: the result of the most recent calculation;
• computer: the computer type;
• eps: the machine epsilon;
• flops: the number of floating point operations;
• i, j: square root of -1;
• inf: "infinity", the result of 1/0, for instance;
• NaN: "Not a Number", the result of 0/0, for instance;
• pi
• realmax, realmin: largest and smallest representable reals;
• version: MATLAB version number.

To see what variable names you've defined, type

who

whos

Variables Are Vectors

We said that MATLAB treats all its variables as though they were vectors (or matrices, or general arrays, but we won't talk about that yet!). So MATLAB considers the variable x we defined earlier to be a vector of length 1. How would we define a longer vector?

The easiest way to define a vector is to list its values inside of square brackets, and separated by spaces or commas:

x = [ 0, 1, 3, 6, 10 ]

x = 1 : 10
x = 50 : 5 : 100
x = linspace ( -1, +1, 25 )
x = linspace ( pi, 2*pi, 11 )

Use the colon notation to set up vectors with whole number entries; use linspace when filling in lots of values between two limits.

Once you've set up your vector, you can change an individual entry, or reference it in a formula, by using index notation. The first element of an array x has index 1 and is named x(1), and so on. Thus, if I want to alter the third element of evens, I could say

evens(3) = 7

Exercise: Define a vector x whose elements are 1, 2, 3, ..., 98, 99, 100. Write a program that computes the sum of the entries of x. You will probably want to use a for loop:

        sum = 0.0
        for i =  start : finish
          statements
        end
      

n * ( n + 1 ) / 2 .

If your variable is a vector, then the length function will tell you how many elements are in it. For instance, after the previous exercise, if we had typed

length ( x )

Vector Operations

Suppose we want a table of integers from 0 to 10, their squares and cubes. We could start with

n = 0 : 10

n2 = n * n

n2 = n .* n

Now let's get back to the problem of vector operations. The multiplication, division and exponentiation operators all have two possible forms, depending on whether you want to operate on the arrays, or on the elements in the arrays. In all these cases, you need to use the PERIOD notation to force elementwise operations. For example, we could also have computed n2 using the exponentiation operator as:

n2 = n .^ 2

Vectors will become very important as you go through the labs. It's important to become comfortable with trying to use them cleverly. It's often possible to perform a task in a few lines by thinking in terms of vectors, where the corresponding elementwise code would be clumsy and tedious.

For instance, let's try to print out our square and cube values. We really, really want them to come out in a nice table form, but we don't want to have to work too hard to do this. What can we do?

If we try this:

          for i = 1 : 11
            n(i), n2(i), n3(i)
          end
      

How about telling MATLAB to group the three numbers as a single object? This gets rid of some garbage:

          for i = 1 : 11
            [ n(i), n2(i), n3(i) ]
          end
      

The very best solution for us is

b = [ n; n2; n3 ];
b'

Polynomials

Naturally, many of the functions we will be looking at will be polynomials. The elementwise operation notation can become pretty cumbersome for these formulas. Let's set up a vector of data points:

x = linspace ( 0.0, 10.0, 101 );

y = x.^3 - 2 * x .* x + x;

Because polynomials come up so often, there's a special command called polyval which evaluates a polynomial given a vector of coefficients, and a vector of evaluation points. For instance, we could have written:

y = polyval ( [ 1, -2, 1, 0 ], x );

c = [ 1, -2, 1, 0 ];
y = polyval ( c, x );

A very useful command for polynomials is the roots command, which returns a vector of the (approximate) roots of a polynomial based on its coefficients:

c = [ 1, -2, 1, 0 ];
x = roots ( c );

y = polyval ( c, x );

Script Files

You will often find that you use the same set of commands repeatedly. These commands might represent an algorithm, or simply a useful little task. Rather than typing the commands again and again, you can type them into a file, name that file, and then run the commands simply by invoking the name of the file. The file should have a name that ends with ".m", as in root.m, euler.m and so on, and these files are called M files.

A simple kind of M file is simply a list of commands, called a script file. You'll like script files if you're a bad typist, or if you need to issue the same set of commands several times. You can even use the editor as needed to modify the commands slightly.

Exercise: using the editor, create a script file called tables.m which contains the lines:

n = [ 0 : 10 ];
n2 = n .* n;
n3 = n .* n2;
b = [ n; n2; n3 ];
b'

Function Files

One special use of an M file is to define new functions, that can be named, and then evaluated just like abs(x) and other familiar functions. To do so, we simply have to write a special kind of M file, with the right name and the right first line, and then calculate the value of the function based on the input.

For instance, we might want to define a function which, for any input x, is equal to sin(x)/x. First we need to pick a name for this function - let's try calling it slinky. Now we need to create a file called slinky.m, and here's what goes inside it:

          function result = slinky ( x )
          %  SLINKY computes sin ( x ) / x.
          result = sin ( x ) / x;
      

There are three lines in the file, and we have a comment about each one of them!

• The first line of the file defines the form of the function; it takes a quantity we'll call x as input, and and returns a quantity we'll call result.
• The % line is a comment. Any comments that occur right after the function line can be automatically printed out, in this case, by typing help slinky.
• The assignment statement has a semicolon at the end. Normally, you don't want to get printout from a function. We probably should worry about what happens at x=0, but it turns out to be a little too much trouble to worry about right now.

Exercise: set up the slinky function. Try typing help slinky and verify that you get the reminder. Make a table of the values of the slinky function at 25 values of x from 0 to pi.

Discussion: Does the value of slinky(pi) bother you? It should be zero, right? For that matter, you should check the value of sin(pi), which is not zero either. In "the real world" it's zero, and in Mathematica it is too. But MATLAB is a numerical calculator, not a symbolic one. In MATLAB, not only is sin(pi) not zero, but pi isn't pi and sin isn't sin.

Graphics

One reason we had to jump into vectors today is that we want to do graphics, one of the nicest features of a language like MATLAB. To make simple graphs, we define vectors x and y of independent and dependent data. Doing a sine graph makes it look easy:

x = linspace ( 0, 10, 100 );
y = sin ( x );
plot ( x, y )

Exercise: Plot the slinky function over the domain of -10 to 10.

c = [ 1, -2, 1, 0 ];
x = linspace ( -2, 4, 50 );
y = polyval ( c, x );
plot ( x, y )

Suppose we also want to plot the x axis, to see where this polynomial is zero. We need the hold on command, which keeps the output of previous plot commands in the plot box, followed by a plot command where we give the vectors x and y explicitly:

hold on
plot ( [ -2, 4 ], [ 0, 0 ] )

Since we're starting a chapter on differential equations, let's try to use graphics to compare the analyze the results of Euler's method. Consider the initial value problem:

y'(x) = 3 * cos(x) * cos(3*x) - sin(x) * sin(3*x)
y(0) = 0.0;

Given that we know the value at x = 0, we will try to estimate the solution value at the 100 points 0.1, 0.2, ..., 10. Given an approximate solution y_n at x_n, Euler's method approximates the solution at x_(n+1) by:

x_(n+1) = x_n + dx
y_(n+1) = y_n + dx * f ( x_n, y_n)

The following exercise might be a good time to try using a script file. Put everything into the script file except for the plotting commands. Do those interactively!

Exercise: Using the sketch below, try to numerically solve the ODE.

dx = ?;
x(1) = ?;
y(1) = ?;
for i = 1 : ?
x(i+1) = x(i) + dx;
y(i+1) = y(i-1) + dx * (derivative);
end

Assignment

Consider the initial value problem

y'(x) = y
y(0) = 1;

        Steps       Computed Y(2.0)    Error
          10        ____________     ____________
         100        ____________     ____________
        1000        ____________     ____________