Introduction To MATLAB Programming

Root-Finding

Now that you are familiar with MATLAB® and its basic functionalities, you will learn how to use MATLAB to find the roots of equations, and specifically, nonlinear equations. We will walk through using Newton’s method for this process, and step through multiple iterations of Newton’s method in order to arrive at a final solution.

Because Newton’s method is an iterative process, you will also learn how to construct two types of logical loops: for loops and while loops. These loops will be useful tools not just for the purpose of using Newton’s method, but also for future use in writing code that can handle more complicated operations. In addition to learning how to construct loops, you will be introduced to plotting in MATLAB and saving code in a file for future or frequent use.

The secant method is presented as an alternative to Newton’s method, and you will be asked to use MATLAB in comparing the convergence of solutions—essentially, the speed at which they can approximate a solution—from these respective methods. When completing this portion of the unit, you will find the need to store multiple values for comparison, and this offers the opportunity to make use of lists and sublists through sub-indexing. Since discrepancies in syntax are important when it comes to using and operating on lists and matrices in MATLAB, you will be asked to evaluate a variety of expressions that help you understand these differences.

Warm-up

Here are some warm-up problems:

  • Create the list of numbers from 3 to 10
  • Create the list of even numbers from –20 to 20
  • Create the decreasing list of numbers from 100 to 90

Remember variables? There’s one variable that is special: The ans variable holds the last calculated value that was not placed into a variable. This is nice when a command you give MATLAB® returns a value that you realize is important, but forgot to assign into a variable:

>> sin(pi/2)
ans =
      1
%%oops! I meant to save the "return value" in the variable x
>> x=ans
x =
      1
>>

Now x holds the answer, 1. Of course you could also re-issue the command with an assignment, but commands can take long to run, and there may be other reasons why you do not want to re-issue the command.

Newton's Method

Many mathematical problems involve solving equations. And while linear equations can be solved rather easily, nonlinear ones cannot. A nonlinear equation can always be written as

 

f(x)=0(1)

 

For a suitably chosen function f. For example, if we want to find x so that tanh(x)=x/3, we could instead chose f(x)=tanh(x)x/3 and solve equation (1). Finding solutions to (1) is called “root-finding” (a “root” being a value of x for which the equation is satisfied).

We almost have all the tools we need to build a basic and powerful root-finding algorithm, Newton’s method*. Newton’s method is an iterativemethod. This means that there is a basic mechanism for taking an approximation to the root, and finding a better one. After enough iterations of this, one is left with an approximation that can be as good as you like (you are also limited by the accuracy of the computation, in the case of MATLAB®, 16 digits).

Iterative methods entail doing the exact same thing over and over again. This is perfect for a computer.

The actual iteration starts from an approximation xn at the nth step and defines the next one, xn+1 :

xn+1=xnf(xn)f(xn)(2)

 

Notice that not only do we need to evaluate the function f at xn, but also we need to evaluate the derivative, f(xn). This could be a problem if the derivative is unknown, or complicated to compute (and there are other methods to use in that case).

For example, we look at a function for which there is no formula for the solution: f(x)=tanh(x)x/3 (we are looking for a non-zero solution). Since tanh(x)=1cosh2(x)=sech2(x), the update rule becomes:

xn+1=xntanh(xn)xn/3sech2(xn)1/3.(3)

 

We can do this manually by starting somewhere (not zero):

>> x=2
x =
     2

and doing one step of the iteration by overwriting the old value with the new one:

>> x=x – (tanh(x) - x/3) / (sech(x)^2 – 1/3)
x =
     3.1320

To see what happens when we do this many times, we would need to tell MATLAB to do it many times. We could type in the line many times, but that would be very annoying and time-consuming. Instead, MATLAB keeps a “Command History” which you can both see (turn it on in the “Desktop” menu) and access by pressing the  and  keys on your keyboard.

Exercise 4. How many iterations do you need until the solution seems to converge? What did it converge to? Tell MATLAB to show the full accuracy and find the 16-digit number that the iteration converges to.

While using the command history is faster than typing in the same command many times, there is an even better way: tell MATLAB to do it many times with either a for loop or a while loop.

Both types of loops tell MATLAB to perform a set of commands over and over again. The idea of a for loop is that you are telling MATLAB to perform a set of command for every element of a list that you supply. For example:

for k=3:10 % for every element in 3:10
     k^2   % calculate _that element_ squared 
end

will iterate over the values 3 to 10 and perform each of the expressions between the for and the end. At each iteration, the variable, k in this case (but can be whatever) will be assigned the appropriate value.

Note 1. When you type code into MATLAB you might make a mistake and want to start over. Press Ctrl C to tell MATLAB to erase the current line and start a new one. Once you start writing a for loop you will not get the command prompt (>>) when you press Enter . This is MATLAB‘s way of telling you that you are in the middle of writing a command. Once you type end and hit Enter , MATLAB will execute the loop and you’ll

get the prompt back. If you started a loop and want to stop in the middle and start over, press 
Ctrl C to get the prompt.

So we can write

% notice how after a % one can write comments which can be used to
% explain the code inline.
% it will actually not bother Matlab if you copy the comments
% together with the code
x=2;                                   % This is our first guess
for k=1:10                             % we will iterate 30 times
     x=x–(tanh(x)–x/3)/(sech(x)^2–1/3) % we do not use k. That is OK.
end

To see the functions we are working with, we can do a nice little plot:

x=linspace(0,pi,100); % a list of 100 points equally spaced between 0
                      % and pi
plot(x,tanh(x),x,x/3) % plot two curves: tanh(x) vs. x and x/3 vs. x

You can read more about the functions we used today using the help command. For example:

>> help tanh
     TANH Hyperbolic tangent.
          TANH(X) is the hyperbolic tangent of the elements of X.

          See also atanh.

          Overloaded methods:
               codistributed/tanh
          Reference page in Help browser
               doc tanh

Notice that in MATLAB help the commands and variables are always written in CAPITAL LETTERS, this is to make them stand out. In practice, the built-in MATLAB commands are all in lower-case. Your own variables (and functions) may use either, and MATLAB is case-sensitive:

>> a=1
a =
     1
>> A
??? Undefined function or variable 'A'.

Exercise 5. Using root-finding calculate ROf course, MATLAB has the function sqrt and also the power function as we saw in the previous lecture. But pretend that it did not. What is RFind a simple function f (that doesn’t use the square-root function) so that f(2)=0(There are several options, so if you don’t manage with one option, try another!) Find 2 like a Babylonian†. How many iterations do you need to get an answer that is 1e–15 from the answer given by MATLAB‡? Note: this problem will require you to use a pencil and paper. You will need to differentiate, divide and simplify a fraction before you type your code in MATLAB.

Notice that the starting point is important. Find starting points that converge to each of ±2.

Exercise 6. As in the previous problem, write an iterative code that will find 1R. You might think that you will have to divide in your code, but that is not true; it depends on the function that you use. Again, if at first you do not succeed, try a different function. Simplify the formula so that it does not need division, and then implement the code to find 1/101. You will need to start close to the answer for the method to converge.

The Secant Method

While Newton’s method is fast, it has a big downside: you need to know the derivative of f in order to use it. In many “real-life” applications, this can be a show-stopper as the functional form of the derivative is not known. A natural way to resolve this would be to estimate the derivative using

 

f(x)f(x+ϵ)f(x)ϵ(1)

 

for ϵ1. The secant method uses the previous iteration to do something similar. It approximates the derivative using the previous approximation. As a result it converges a little slower (than Newton’s method) to the solution:

 

xn+1=xnf(xn)xnxn1f(xn)f(xn1).(2)

 

Since we need to remember both the current approximation and the previous one, we can no longer have such a simple code as that for Newton’s method. At this point you are probably asking yourself why we are not saving our code into a file, and it is exactly what we will now learn how to do.

Coding in a File

Instead of writing all your commands at the command prompt, you can type a list of commands in a file, save it and then have MATLAB® “execute” all of the commands as if you had typed them into the command prompt. This is useful when you have more than very few lines to write because inevitably you are bound to make a small mistake every time you write more than 5 lines of code. By putting the commands in a file you can correct your mistakes without introducing new ones (hopefully). It also makes it possible to “debug” your code, something we will learn later.

For guided practice and further exploration of how to debug, watch Video Lecture 6: Debugging.

MATLAB files have names that end with .m, and the name itself must comprise only letters and numbers with no spaces. The first character must be a letter, not a number. Open a new file by clicking on the white new-file icon in the top left of the window, or select from the menu FileNewScript. Copy the Newton method code for tanh(x)=x/3 into it. Save it and give it a name (NewtonTanh.m for example). Now on the command prompt you “run” the file by typing the name (without the .m) and pressing Enter .

A few points to look out for:

  • You can store your files wherever you want, but they have to be in MATLAB’s “search path” (or in the current directory). To add the directory you want to the path select FileSet path… select “Add Folder”, select the folder you want, click “OK” then “Save”. To check if your file is in the path you can type which NewtonTanh and the result should be the path to your file.
  • If you choose a file-name that is already the name of a MATLAB command, you will effectively “hide” that command as MATLAB will use your file instead. Thus, before using a nice name like sum, or find, or exp, check, use which to see if it already defined.
  • The same warning (as the previous item) applies to variable names, a variable will “hide” any file or command with the same name.
  • If you get strange errors when you try to run your file, make sure that there are no spaces or other non-letters in your filename, and that the file is in the path.
  • Remember that after you make changes to your file, you need to save it so that MATLAB will be aware of the changes you made.

For guided practice and further exploration of how to use MATLAB files, watch Video Lecture 3: Using Files.

Exercise 7. Save the file as SecantTanh.m and modify the code so that it implements the Secant Method. You should increase the number of iterations because the Secant Method doesn’t converge as quickly as Newton’s method.

 

Notice that here it is not enough to use x like in the Newton’s method, since you also need to remember the previous approximation xn1Hint: Use another variable (perhaps called PrevX).

Convergence

Different root-finding algorithms are compared by the speed at which the approximate solution converges (i.e., gets closer) to the true solution. An iterative method xn+1=g(xn) is defined as having pth order convergence if for a sequence xn where limnxn=α exists then

 

limn|xn+1α||xnα|p=L0.(3)

 

Newton’s method has (generally) second-order convergence, so in Eq. (3) we would have p=2, but it converges so quickly that it can be difficult to see the convergence (there are not enough terms in the sequence). The secant method has a order of convergence between 1 and 2. To discover it we need to modify the code so that it remembers all the approximations.

The following code, is Newton’s method but it remembers all the iterations in the list x. We use x(1) for x1 and similarly x(n) for xn:

x(1)=2;   % This is our first guess, put into the first element of x
for n=1:5 % we will iterate 5 times using n to indicate the current
          % valid approximation
     x(n+1)=x(n)-(tanh(x(n))-x(n)/3)/(sech(x(n))^2-1/3);  %here we
                                % calculate the next approximation and
                                % put the result into the next position
                                % in x.
end
x % sole purpose of this line is to show the values in x.

The semicolon (;) at the end of line 4 tells MATLAB not to display the value of x after the assignment (also in line 1. Without the lonely x on line 9 the code would calculate x, but not show us anything.

After running this code, x holds the 6 approximations (including our initial guess) with the last one being the most accurate approximation we have:

x =
     2.0000    3.1320    2.9853    2.9847    2.9847    2.9847

Notice that there is a small but non-zero distance between x(5) and x(6):

>> x(6)-x(5)
ans =
     4.4409e-16

This distance is as small as we can hope it to be in this case.

We can try to verify that we have second order convergence by calculating the sequence defined in Eq. (3). To do that we need to learn more about different options for accessing the elements of a list like x. We have already seen how to access a specific element; for example to access the 3rd element we write x(3). MATLAB can access a sublist by giving it a list of indexes instead of a single number:

>> x([1 2 3])
ans =
     2.0000    3.1320    2.9853

We can use the colon notation here too:

x(2:4)
ans =
     3.1320    2.9853    2.9847

Another thing we can do is perform element-wise operations on all the items in the list at once. In the lines of code below, the commands preceding the plot command are executed to help you understand how the plot is generated:

>> x(1:3).^2
ans =
     4.0000    9.8095    8.9118
>> x(1:3)*2
ans =
     4.0000    6.2640    5.9705
>> x(1:3)-x(6)
ans =
    -0.9847    0.1473    0.0006
>> x(2:4)./(x(1:3).^2)
ans =
     0.4002    0.0018    0.0387
>> plot(log(abs(x(1:end-2)-x(end))),log(abs(x(2:end-1)-x(end))),'.'))

The last line makes the following plot (except for the green line, which is y=2x):

Graph of a line and roots to verify second-order convergence of Newton's method.

MATLAB can calculate roots through Newton’s method, and verification of convergence is graphed.

The main point here is that the points are more or less on the line y=2x, which makes sense:

Taking the logarithm of the sequence in (3) leads to

 

log|xn+1α|logL+plog|xnα|(4)

 

for n1, which means that the points (log|xnα|,log|xn+1α|) will converge to a line with slope p.

The periods in front of */, and ^ are needed (as in the code above) when the operation can have a linear algebra connotation, but what is requested is an element-by-element operation. Since matrices can be multiplied and divided by each other in a way that is not element-by-element, we use the point-wise version of them when we are not interested in the linear algebra operation.

Exercise 8. Internalize the differences between the point-wise and regular versions of the operators by examining the results of the following expressions that use the variables A=[1 2; 3 4]B=[1 0; 0 2], and C=[3;4]. Note: some commands may result in an error message. Understand what the error is and why it was given.

  • A*B vs. A.*B vs. B.*A vs. B*A
  • 2*A vs. 2.*A
  • A^2 vs. A*A vs. A.*A vs. A.^2 vs. 2.^A vs. A^A vs. 2^AThe last one here might be difficult to understand…it is matrix exponentiation.
  • A/B vs. A\B vs. A./B vs. A.\B
  • A*C vs. A*C' vs. C*A vs. C'*A
  • A\C vs. A\C' vs. C/A vs. C'/A

Homework 2. Modify your secant method code so that it remembers the iterations (perhaps save it in a new file?). Now plot the points that, according to (4) should be on a line with slope p. What is p?

More Sub-Indexing

We have seen how one can access a subset of a list by providing a list of desired positions:

>>  A=rand(1,5)
A =
     0.6430    0.5461    0.5027    0.0478    0.2289
>> A([2 3 1 1])
ans =
     0.5461    0.5027    0.6430    0.6430

There are a few more extensions of this:

It can be used to modify part of a matrix:

%with A as before: 
>> A([1 2])=3+A([3 4])
A =
     3.5027    3.0478    0.5027    0.0478    0.2289
%or even:
>> A([3 4])=1
A =
     3.5027    3.0478    1.0000    1.0000    0.2289 

Additionally, this works for matrices and submatrices as well:

>> A=magic(4)
A =
     16      2     3     13
      5     11    10      8
      9      7     6     12
      4     14    15      1
>> A([2 3], [1 4])
ans =
     5     8
     9    12
>> A(1:2,3:4)=1
A =
     16     2     1     1
      5    11     1     1
      9     7     6    12
      4    14    15     1

The keyword end will evaluate to the size of the dimension of the matrix in which it is located:

>> A=magic(4);
>> A(end,end)
ans =
     1
>> A([1 end],[1 end])
ans =
     16    13
      4     1
>> A([1 end/2], [2, end-1])
ans =
      2     3
     11    10
>> A(2,1:end)
ans =
     5    11    10     8
>> A(1:end,3)
ans =
      3
     10
      6
     15
%  1:end is so useful that it has an even shorter notation, :
>> A(:,1)
ans =
     16
      5
      9
      4
>> A(4,:)
ans =
     4    14    15     1

Finally, a matrix can be accessed with a single index (as opposed to with two) and this implies a specific ordering of the elements (rows first, then columns):

>> A(5)
ans =
     2
>> A(4:10)
ans =
     4     2    11     7    14     3    10
>> A(:)
ans =
     16
      5
      4
      2
     11
      7
     14
      3
     10
      6
     15
     13
      8
     12
      1

Exercise 9: Practice some of these methods for accessing matrix elements with the following exercises:

    • Create a matrix of size N×N that has ones in the border and zeros inside. For example if N=3 the matrix can be created with
>> A=ones(3,3); A(2,2)=0
A =
     1     1     1
     1     0     1
     1     1     1

Make this construction depend on N and work for any positive integer N2

  • Create a 5×5 matrix whose rows are (1:5)
  • Extract the diagonal of a given matrix without using diag (you may use size)
  • Flip a given matrix horiztonally. Vertically? Do not use fliplr or flipud
  • Extract the anti-diagonal of a given matrix
  • Extract the anti-diagonal, without first flipping it (Hint: use single index access)