MATH/APPM 4650
Class notes, January 23

Roundoff and Introduction to Mathematical Software




More on roundoff error

The same sort of thing happens with the quadratic equation. a x² + bx + c = 0.

x = [ -b +- sqrt(b² - 4ac) ] / [2a].
If b is large and positive, then the term with the plus sign involves subtracting two nearly equal numbers.
(The term with the minus sign is fine.)

Example: four significant digits. x² + 62.10x + 1 = 0.

x = [-62.10 +- sqrt(62.10² - 4)] / 2.
Exact solutions (to four places) are -0.1611 and -62.08.

62.10² = 3856.41 which rounds to 3856.
So 62.10² - 4 = 3852.
sqrt(3852) = 62.06448258, which rounds to 62.06.
(62.06-62.10)/2 = -0.2000.     (-62.06-62.10)/2 = -62.08.
First root is WAY off, second root is correct to tolerance.


Rationalize: x = [-2c] / [-+ sqrt(b²-4ac) + b].
Now compute sqrt(62.10²-4) = 62.06, so
-2/(62.06+62.10) = -0.1611     -2/(-62.06+62.10) = -50.00.
First root is correct to tolerance, second root is WAY off.


Hence neither form of the solution gives a good algorithm for solving the quadratic equation.
We need to use conditioning depending on the data to get both roots.

Lesson: Exact formulas can sometimes be worse than approximate techniques!
Lesson: Mathematically equivalent formulas are not numerically equivalent!

---

In general, roundoff error is unpredictable.
Even a few additions and multiplications can accumulate roundoff error.

Therefore one should always minimize the number of computations one does.
For example,
a ( b + c ) = ab + ac mathematically.
However, the left side involves two computations, while the right side involves three.
Roundoff error is more likely to accumulate in the right side.
So we should use the left side for numerics.

When evaluating a polynomial, the same things arise.
Instead of computing ax² + bx+c directly (which involves five operations),
we should compute x(ax+b) + c (which involves four operations).
The difference seems minor, but it adds up.

See the results of the Maple worksheet.

---

General rules for roundoff:

• Adding and multiplying will produce small roundoff errors.
• Subtracting numbers that are nearly equal gives HUGE roundoff errors.
• Dividing by small numbers magnifies errors GREATLY.
• Roundoff error is hard to predict. Always assume the worst.
• The more computations, the more roundoff errors. The fewer the better!
• Order of operations always matters numerically. Equivalent formulas often produce different numerical results.

---

Mathematical Software and programming overview

These links will help you get started with basic computations and commands in the three software programs.
• A brief Maple tutorial by John Pais: http://kinetigram.com/micds/Maple10Tutorial/quickintroSum05.mws (Open Maple, then select "Open URL" from the File menu and copy & paste this link into it. Click the arrows to open each section up. More tutorials available at the official site.)
• Official MATLAB tutorial: http://www.mathworks.com/academia/student_center/tutorials/prerequisites.html (Follow along with suggested commands to get a feel for the format.) There are also video tutorials and interactive tutorials at the official site.
• Official Mathematica tutorial: http://library.wolfram.com/conferences/devconf99/withoff/index2.html (Copy commands into your window and follow along. More help available at the official site.)

All the programs also have extensive help from within the software.


The university also maintains a list of common commands in the three major programs..


For the most part, we will want to write programs.

The programs we will write will have the following basic structure:

1. Assign initial variables.
2. Begin an iteration.
3. Modify the variables in some way.
4. Increment an index and/or check an error.
5. Check a condition, and either end the iteration or continue.
6. When the iteration ends, output one or more final variables.

Example: A pseudocode program to approximate Euler's constant e = 2.71828...

We use the clever method: e = 1 + 1 + 1/2 * (1 + 1/3 * (1 + 1/4 * (1 + 1/5 * (1 + 1/6 * (... + 1/29 * (1 + 1/30) ... ) ) ) )
EULER = 1
N = 3
FOR J FROM 1 TO N DO
      EULER = EULER * 1/(N-J+1)
      EULER = EULER + 1
END DO
PRINT (EULER)

Let's see what this program is doing.

1. Set the initial approximation to 1.
2. Set the number of computations to do: N. (Observe that we can change N easily to run the same program for 3 steps, 30 steps, or 300 steps.) Let's see what happens when N = 3.
3. Start the iteration. This will run N steps. The variable J is a counter. It will run the indented instructions with J=1, then with J=2, then get out of the loop at "END DO"
4. First time through: N=3, J=1 and EULER=1. Compute 1/(N-J+1) which is 1/3. So EULER is now 1/3.
5. Add 1 to EULER. So EULER is now 1+1/3.
6. Go back, increment J. If J is less than or equal to N, keep going. (It is.)
7. Second time through: N=3, J=2, and EULER = 1+1/3. Compute 1/(N-J+1) which is 2. So EULER is now 1/2*(1+1/3).
8. Add 1 to EULER. So EULER is now 1+1/2*(1+1/3).
9. Go back, increment J. If J is less than or equal to N, keep going. (It is.)
10. Third time through: N=3, J=3, and EULER = 1+1/2*(1+1/3). Compute 1/(N-J+1) which is 1. So EULER is now 1*(1+1/2*(1+1/3)).
11. Add 1 to EULER. So EULER is now 1+1*(1+1/2*(1+1/3)).
12. Go back, increment J. If J is less than or equal to N, keep going. It's not! So get out of the loop and execute the next command.
13. This prints the result.
14. A hand computation verifies that the answer should be 8/3 = 2.666..., so we check that the program gives us this.
15. If so, we can trust the output, and run it for N=10.
16. The answer we get is 2.718281801, which is correct to eight places.

Maple code for this program. (Observe the use of "evalf" to get floating point numbers, and the colons to suppress output at each step. To show output, use semicolons instead. Exactly the opposite of the others!)
Enter each line as written, and press SHIFT+ENTER after each line. Then press ENTER at the end to run it.
> Euler:=1:
  N:=3:
  for J from 1 to N do
     Euler := evalf(Euler * 1/(N-J+1)):
     Euler := evalf(Euler + 1):
  end do:
  print(Euler):


MATLAB code for this program. (Observe the semicolons to suppress output at each step. To show output, remove semicolons.)
Enter each line as written, and press SHIFT+ENTER after each line. Then press ENTER at the end to run it.
Alternatively, you can copy the entire thing into an Editor window and press F5 to run it. (This makes it easier to change values.)
  Euler = 1;
  N = 3;
  for J=1:N
      Euler = Euler*1/(N-J+1);
      Euler = Euler + 1;
  end
  Euler


Mathematica code for this program. (Observe the use of N[ ... ] to get floating point numbers out and semicolons to suppress output. To show output, remove semicolons.)
Enter each line as written, and press ENTER after each line. Then press SHIFT+ENTER at the end to run it. (Exactly the opposite of the others!)
  Euler=1;
  M=3;
  For[J=1,J<=M,J++,
       {Euler=N[Euler*1/(M-J+1)],
       Euler=N[Euler+1]}
    ]
  Print[Euler];



Example: A pseudocode program to estimate sqrt(2) = 1.41421...

We use the popular algorithm x_new = x_old/2 + 1/x_old. The idea is to keep applying this algorithm until we get within the desired tolerance, or we use too many iterations (whichever comes first.)

X = 1
MAXITER = 2
COUNT = 1
TOL = 0.001
ERROR = TOL + 1
WHILE ( (ERROR > TOL) AND (COUNT <= MAXITER) ) DO
   X = X/2 + 1/X
   ERROR = ABS(X*X-2)
   COUNT = COUNT + 1
END DO
IF (ERROR > TOL) THEN
   PRINT("Failed to reach tolerance")
ELSE
   PRINT("Reached tolerance")
END IF
PRINT(X)


What is this program doing?

1. Start with an initial guess x = 1.
2. Set MAXITER = 2. No matter what, we will not do more than two computations.
3. Set COUNT = 1. This will keep track of how many computations we did.
4. Set TOL = 0.001. We'd like the square of our final answer to be within 0.001 of 2.
5. Set ERROR = TOL + 1. Set the initial error to be anything greater than TOL.
6. Now enter the WHILE loop. If both conditions are satisfied, run the computations inside.
7. Is ERROR larger than TOL? Yes. Is COUNT less than or equal to MAXITER? Yes.
8. First time through: So the new x is the average of the old x and 2/x. So X is 1.5 now.
9. Compute the difference between x² and 2. That's our new ERROR. It's 0.25 now.
10. Increase COUNT. It's 2 now.
11. Now go back and check the conditions again.
12. Is the ERROR larger than TOL? Yes. Is COUNT less than or equal to MAXITER? Yes.
13. Second time through: The new X is 1.5/2 + 1/1.5 = 1.41666...
14. The new ERROR is 0.007.
15. The new COUNT is 3.
16. Go back and check the conditions again.
17. Is the ERROR larger than TOL? Yes. Is COUNT less than or equal to MAXITER? No!
18. Since both conditions are not satisfied, we do not go through the loop again.
19. Now check why the loop stopped. If ERROR is larger than TOL, the loop stopped because COUNT reached MAXITER.
20. Print a bad message.
21. If ERROR had been smaller than or equal to TOL, we would have printed a good message.
22. Finally, print the answer.

Maple code for this program:
> x:=1:
  Maxiter:=2:
  Count:=1:
  Tol:=0.001:
  Err:=Tol+1:
  while ((Err>Tol) and (Count<=Maxiter)) do
     x:=evalf(x/2+1/x):
     Err:=evalf(abs(x^2-2)):
     Count:=Count+1;
  end do:
  if (Err>Tol) then
     print(`Failed to reach tolerance`)
  else
     print(`Reached tolerance`)
  end if:
  print(x):


MATLAB code for this program:
  x=1;
  Maxiter=2;
  Count=1;
  Tol=0.001;
  Err=Tol+1;
  while (Err>Tol) && (Count<=Maxiter)
     x=x/2+1/x;
     Err=abs(x^2-2);
     Count=Count+1;
  end
  if Err>Tol
     disp('Failed to reach tolerance');
  else
     disp('Reached tolerance');
  end
  x


Mathematica code for this program:
  x=1;
  Maxiter=2;
  Counter=1;
  Tol=0.001;
  Err=Tol+1;
  While[ ( (Err> Tol) && (Counter<=Maxiter) ),
    {x=N[x/2+1/x];
     Err=Abs[x^2-2];
     Counter += 1}
  ]
  If[Err>Tol,
     Print["Failed to reach tolerance"],
     Print["Reached tolerance"]
  ]
  Print[x]

---

This should get you through most of the basic programming commands.

All of the programs in the textbook, for each of the formats, are on the author's web site.

• Maple versions. The simplest thing to do is to download "Maple-Programs.zip". Extract the files somewhere, and use the ones from the folder "mws-6". Open "Classic Worksheet Maple" to use these (a regular Maple session will choke on them). Execute the program by repeatedly pressing "Enter" and change the values as desired.
• MATLAB versions. Download "Matlab-Programs.zip". Open the programs in an Editor window and press F5 to run them. You will need the Symbolic Math Toolbox (for some reason).
• Mathematica versions. Again, you can just download "Mathematica-Programs.zip" and extract the files somewhere. Execute the programs by pressing "Shift-Enter" and enter values as desired. Note that some of these programs will not end properly. You may need to force a quit!