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
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;Is this how you set up a formula F(X)?
y = 4;
z = sqrt ( x^2 + y^2 )
x = 1;Apparently not, since f doesn't change when x does. We'll see how to set up a formula f(x) later.
f = 1 / ( 1 + x^2 )
x = 4;
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 = 1but, first of all, there's really no reason to store each result in a new variable, and secondly, we can use the do command to repeat this for us as often as we like:
x2 = 1 + 1 / x1
x3 = 1 + 1 / x2
x4 = 1 + 1 / x3
x = 1(Do you know what the limiting value of x will be?)
for i = 1 : 10
x = 1 + 1 / x
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.
We never thought about what goes on when we say
x = 1to MATLAB. In a programming language, you would probably have to declare that x was a variable, and that it had a particular arithmetic type (say integer or real or complex) and that it was a scalar (just one number), or a vector (a string of numbers), or an array (a table). We did none of that.
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 epsilonWe weren't careful. The value of epsilon that we just computed is incorrect. Can you tell me why? Can you fix the code, or fix the value of epsilon?
You can test your value of epsilon in two ways:
( 1 + epsilon ) - 1
One of the simplest commands in MATLAB has the form:
yUsually, this command returns the value of the variable, or built-in function, or user-defined formula. So the first idea to keep in mind is this:
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:
To see what variable names you've defined, type
whofor a short list, or
whosfor a lot of information.
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 ]This is not practical if you want to set up a whole lot of data. In that case, assuming your values are regularly spaced, you can use one of the following commands:
x = 1 : 10The first command assumes an increment of 1 between the values. In the second command, an increment of 5 is specified. In the third command, 25 equally spaced values are generated, while the fourth generates 11.
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 endCheck your result, given that the sum of the integers from 1 to n is:
n * ( n + 1 ) / 2 .Now use the start:increment:finish notation, and modify your program slightly, to compute the sum of the odd elements in the array (that is, the elements in entries 1, 3, 5, and so on).
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 )we should get back the value 100. This works whether x is a row vector or column vector. If you don't even remember whether your vector is a row or column vector, you can issue the size command, which will list the number of rows and columns.
Suppose we want a table of integers from 0 to 10, their squares and cubes. We could start with
n = 0 : 10but now we'll get an error when we try to multiply the entries of n by themselves.
n2 = n * nRealize that MATLAB thinks you're dealing with vectors, and the default multiplication operation with vectors is vector multiplication. In order to do element-by-element multiplication, we need to place a period in front of the operator:
n2 = n .* nNow we can define n3 in a similar way.
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 .^ 2These problems never come up with addition or subtraction; nor do they occur with division or multiplication by a scalar.
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) endwe get all kinds of MATLAB garbage that makes the table hard to read.
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 ];The first command stores each vector as a row in b. The second command transposes b, and that's what we want to see!
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 );and now consider the problem of evaluating the polynomial x^3-2*x^2+x at those points. We might type
y = x.^3 - 2 * x .* x + x;We have to be careful to put in periods in both the exponentations and the multiplications. Notice that where I multiply by 2, I don't have to use the elementwise operator!
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 );Or it might even make sense to save the polynomial coefficients in a vector first:
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 ];Following these two commands, what would be the result of:
x = roots ( c );
y = polyval ( c, x );
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 ];Now type the command tables and see what happens.
n2 = n .* n;
n3 = n .* n2;
b = [ n; n2; n3 ];
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!
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.
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.Now, let's try to plot the polynomial, x^3-2*x^2+x:
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 onOf course, the hold off command will cause the next plot command to start with a fresh screen.
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)whose solution is y(x) = cos(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 yn at xn, Euler's method approximates the solution at x(n+1) by:
x(n+1) = xn + dx
y(n+1) = yn + dx * f ( xn, yn)
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 = ?;Using the same set of x values, define the vector y2 to be the value of the true solution at those points. Using the hold on command, compare the true and computed solution on one plot. Satisfied?
x(1) = ?;
y(1) = ?;
for i = 1 : ?
x(i+1) = x(i) + dx;
y(i+1) = y(i-1) + dx * (derivative);
Consider the initial value problem
y'(x) = yUse Euler's method to approximate the solution at x = 2. Do this three times, using 10, 100 and 1000 steps. The exact solution at the end point is exp(2.0). Compute the error in each of your approximate solutions, that is, the difference between your computed value of Y(2.0) and the exact value. Record your results in the following table:
y(0) = 1;
Steps Computed Y(2.0) Error 10 ____________ ____________ 100 ____________ ____________ 1000 ____________ ____________Make some intelligent one-sentence comment about the behavior of the error. Mail the table and your comment to me.