1
Math Homework 471, Fall 2006 Assigned: Friday, September 7, 2007 Due: Friday, September 14, 2007 Section 002
This document is intended to show what a good solution to the first homework assignment might look like. Note in particular the following points of style: There is a cover page with my name and section number. I have written in complete, English sentences. I boxed my answers whenever it was appropriate. (I admit that the boxes are not typeset very well, but we all have our flaws.) The answers are both complete and concise. I included the Matlab code for my program and a transcript of the output. The Matlab code has comments. The tables are neatly printed with the required number of rows, decimal places, etc. The problems are all labeled, and the pages are all in order.
1
2
(1) (Finite precision numbers) The floating point representation of a real number takes the form x = (0.a-1 a-2 . . . a-n ) e , where a-1 = 0, -M e M . Suppose that = 2, n = 4, and M = 5. (a) Find the smallest and largest positive numbers that can be represented in this floating point system. Give your answers in decimal form. The largest positive number that can be represented in this floating point system is +(0.1111)2 25 = 30. The smallest positive number is +(0.1000)2 2-5 = 0.015625 (Recall that the first digit to the right of the decimal must be nonzero.) (b) Find the floating point number in this system that is closest to The number closest to 2 = 1.4142 is 2.
+(0.1011)2 21 = 1.375.
3
(2) (Programming) Implement the following algorithm for converting a decimal number to its binary equivalent x = (an an-1 . . . a0 )2 . Step 1: x0 := x; j := 0 Step 2: while xj = 0, do aj := remainder of integer division xj /2 xj+1 := quotient of integer division xj /2 j := j + 1 end while The Matlab functions rem, mod, and floor may be useful. Make sure your algorithm prints the result out in the correct order. Apply the algorithm to convert the numbers x = 34 and x = 2006 into binary. The following Matlab function accepts a natural number x and prints its binary expansion. I've implemented some error checking and added commentary since it's good programming practice. function PrintBinary(x) if (x <= 0), error( 'The input was not positive.' ); end if (x ~= floor(x)), error( 'The input was not an integer.' ); end j = 1; % Set the position counter to one % since Matlab indexes arrays starting % from one (not zero, like C)
while (x > 0) a(j) = mod(x, 2); x = floor(x / 2); j = j + 1; end a = a(:)'; a = fliplr( a ); a % Force a to be a row vector % Reverse the entries so they print % in the correct order % Print out the answer % Find and store the LSB of x % Remove the LSB from x
Here is a transcript of the program output, edited slightly for legibility. >> PrintBinary( 34 ) 1 0 0 0 1 0 >> PrintBinary( 2006 ) 1 1 1 1 1 0 1 0 1 1 0
4
(3) (Rounding arithmetic) Use three-digit, decimal rounding arithmetic (i.e., = 10 and n = 3) to compute the following sums. Add the numbers by hand in the specified order.
6
(a)
k=1
1 3k
6
(b)
k=1
1 37-k
The first sum equals 0.498 . 0.333 0.111 0.444 0.0370 0.481 0.0123 0.493 0.00412 0.497 0.00137 0.498 The second sum equals 0.499 . 0.00137 0.00412 0.00549 0.0123 0.0178 0.0370 0.0548 0.111 0.166 0.333 0.499
+ = + = + = + = + =
+ = + = + = + = + =
The correct answer 0.499314 = 0.499 in three-digit, decimal, rounding arithmetic. The reason that the first sum is incorrect is that we suffer a loss of precision by adding the small numbers last. In general, it is better to add small numbers first.
5
(4) (Taylor series) (a) Calculate the first three terms of the Taylor series for f (x) = ln x about x = 2. Express the coefficients analytically (not as decimals). The second-degree Taylor polynomial for f (x) = ln x expanded about x = 2 is P2 (x) = ln 2 + 1 (x - 2) - 2
1 2!4 (x
- 2)2 .
(b) The first-order Taylor expansion of f (x) = exp x about x = 0 with exact remainder is f () 2 x 2! for some point between zero and x. Find an analytic expression for . Use this expression to determine the numerical value of when x = 0.5 and x = 0.1. f (x) = 1 + x + The first-order Taylor expansion of the exponential function about zero is e 2 x . 2! Solving this equation for , we obtain an explicit formula: ex = 1 + x + = ln 2(ex - 1 - x)x-2 . When we apply this formula with x = 0.5 we see that = 0.17376 which is clearly between x0 = 0 and x = 0.5. When we apply this formula with x = 0.1, we see that = 0.03361 which is between x0 = 0 and x = 0.1.
6
(5) (Cancellation Errors, p52 #12) Near certain values of x, the following functions cannot be accurately computed using the given formula on account of arithmetic cancellations. Identify the values of x where cancellation occurs (e.g., near x = 0 or when x is large and positive). Propose a reformulation that removes the problem (e.g., using Taylor series, rationalization, trigonometric identities, etc.). (a) f (x) = 1 + cos x (c) f (x) = ln x - ln(1/x) (e) f (x) = 1 - 2 sin2 x (g) f (x) = x - sin x (b) f (x) = + sin x 1 e-x - 2 + 1 - x2 + 4 (d) f (x) = x (f) f (x) = ln(x + x2 + 1) (h) f (x) = ln x - 1
(a) The function f (x) = 1 + cos x has cancelation near x = k where k = 1, 3, . . . . Using a trigonometric identity, we rewrite f (x) = 2 cos2 x.
(b) The function f (x) = e-x + sin x - 1 exhibits cancelation near x = 0 . Simplifying its fourth-order Taylor expansion about zero, we obtain f (x) = 1 - x + = x2 x3 x3 x3 - + + O(x5 ) + x - + O(x5 ) - 1 2! 3! 4! 3!
x2 x3 - + O(x4 ). 2 3
Therefore, for x 0, we have f (x) x2 /2 - x3 /3 . (c) The function f (x) = ln x - ln(1/x) exhibits cancelation near x = 1 . Using properties of logarithms, express f (x) = ln x2 . (d) The function f (x) = x2 + 1 - x2 + 4 shows cancelation when |x| . The following trick solves the problem: f (x) = =
x2 + 1 -
x2 + 4 x2 + 1 + x2 + 4 x2 + 1 + x2 + 4
-3 . x2 +1+ x2 +4
(e) The function f (x) = 1 - 2 sin2 x has problems when x = /4 (and other values). Using a trigonometric identity, f (x) = 1 - 2 1 - cos(2x) = cos(2x) . 2
7
(f) The function f (x) = ln(x + x2 + 1) is difficult to evaluate when x - . Using rationalization and properties of logarithms, x + x2 + 1 x - x2 + 1 f (x) = ln x - x2 + 1 = - ln x - x2 + 1 .
(g) The function f (x) = x - sin x exhibits cancelation near x = 0 . The fifth-order Taylor expansion gives f (x) = x - x - x3 x5 + + O(x7 ) 3! 5! = x3 x5 - + O(x7 ). 3! 5!
For x 0, we have f (x) x3 /3! - x5 (h) /5!. The function f (x) = ln x-1 has problems near x = e . Using properties of logarithms, we write f (x) = ln x . e
8
(6) (Rate of convergence of a sequence, p27 #1) Compute the following limits and determine the corresponding rates of convergence. (a) limn
n-1 n3 +2
Since the numerator is a lower-degree polynomial than the denominator, the limit is zero . Since n-1 n 1 3 = 2, 3+2 n n n and the rate of convergence is O(n-2 ). (b) limn ( n + 1 - n) We may use rationalization to see that 1 1 n+1- n= . 2 n n+1+ n Therefore, the limit is zero , and the rate of convergence is O(n-1/2 ). (c) limn
sin n n
Since | sin n| 1, the limit of the sequence is zero , and the rate of convergence is O(n-1 ).
9
(7) (Rate of convergence of a function) Let f be a function. The following function approximates the derivative of f at a point x when h is small and positive. f (x + h) - f (x) g(x; h) = h Let f (x) = sin x and fix x = 1. Define e(h) = |g(1; h) - cos 1|. For h = 2-n with n = 1, 2, . . . , 8, complete the following table. (Use eight or more decimal places.) n h g(1; h) e(h) e(h)/h e(h)/h2 e(h)/h3 Use this table to decide the rate of convergence of the sequence. Explain your reasoning. Repeat this process for the function f (x + h) - f (x - h) g(x; h) = , 2h which also approximates the derivative of f . The results for the first approximation follow. We see that the g(h) approximates cos 1 = 0.54030231 as h approaches zero since the error e(h) tends to zero. The rate of convergence is O(h) since the fifth column of the table is roughly constant. In other words, the error is proportional to h as h tends to zero.
====================================================================================== n h g(h) e(h) e(h)/h e(h)/h^2 e(h)/h^3 ====================================================================================== 1 0.5 0.312048004 0.22825430 0.45650860 0.91301721 1.82603442 2 0.25 0.430054538 0.11024777 0.44099107 1.76396428 7.05585713 3 0.125 0.486372874 0.05392943 0.43143545 3.45148362 27.61186895 4 0.0625 0.513663206 0.02663910 0.42622560 6.81960963 109.11375410 5 0.03125 0.527067456 0.01323485 0.42351519 13.55248611 433.67955567 6 0.015625 0.533706463 0.00659584 0.42213395 27.01657297 1729.06067012 7 0.0078125 0.537009830 0.00329248 0.42143687 53.94391922 6904.82166034 8 0.00390625 0.538657436 0.00164487 0.42108672 107.79819941 27596.33904960 ======================================================================================
The results for the second approximation follow. As before, we see that the g(h) approximates cos 1 = 0.54030231 as h approaches zero since the error e(h) tends to zero. Now, the rate of convergence is O(h2 ) since the sixth column of the table is roughly constant. In other words, the error is proportional to h2 as h tends to zero. This suggests that the second rule (i.e., central differencing), might be a better way to approximate derivatives than the first rule (i.e., forward differencing).
================================================================================= n h g(h) e(h) e(h)/h e(h)/h^2 e(h)/h^3 ================================================================================= 1 0.5 0.51806945 0.02223286 0.04446572 0.08893143 0.17786286 2 0.25 0.53469172 0.00561059 0.02244235 0.08976940 0.35907758 3 0.125 0.53889637 0.00140594 0.01124751 0.08998006 0.71984047 4 0.0625 0.53995062 0.00035169 0.00562705 0.09003280 1.44052477 5 0.03125 0.54021437 0.00008794 0.00281394 0.09004599 2.88147160 6 0.015625 0.54028032 0.00002198 0.00140702 0.09004929 5.76315424 7 0.0078125 0.54029681 0.00000550 0.00070352 0.09005011 11.52641402 8 0.00390625 0.54030093 0.00000137 0.00035176 0.09005032 23.05288079 =================================================================================
10
(8) (Order of convergence of a sequence) The following sequence converges to 5. Let p0 = 1. For n 0, define 5 1 pn + pn+1 = 2 pn Define the absolute error en = |pn - 5|. Make a table with the following columns, letting n = 0, 1, . . . , 6. (Use eight or more decimal places.)
3 n pn en en /en-1 en /e2 n-1 en /en-1
Use this table to decide the order of convergence of the sequence. Explain your reasoning. Repeat this process for the sequence defined by p0 = 1 which also converges to and pn+1 = p3 + 5 3pn n 3p2 + 5 n for n 0,
5. Consider n = 0, 1, . . . , 4.
The results for the first sequence follow. We see that the iterates do converge to 5 = 2.236069775 and the error converges to zero. The sequence probably converges quadratically (order = 2) , since the fifth column appears to approach a constant. It is clear at least that the sequence does not converge linearly (order = 1) nor cubically (order = 3) since neither the fourth nor sixth column is anywhere near constant. Once the error reaches zero, we obtain the quantity NaN (i.e., Not a Number) on account of a zero divide.
=========================================================================== n p(n) e(n) e(n)/e(n-1) e(n)/e(n-1)^2 e(n)/e(n-1)^3 =========================================================================== 0 1.00000000 1.23606798 1 3.00000000 0.76393202 0.618033989 0.50000000 0.40450850 2 2.33333333 0.09726536 0.127322004 0.16666667 0.21816950 3 2.23809524 0.00202726 0.020842576 0.21428571 2.20310420 4 2.23606890 0.00000092 0.000452899 0.22340426 110.20006787 5 2.23606798 0.00000000 0.000000205 0.22336488 243278.80918904 6 2.23606798 0.00000000 NaN NaN NaN ===========================================================================
The results for the second sequence follow. As before, the iterates converge to 5, and the error converges to zero. The sequence appears to converge cubically (order = 3) since the last column seems to be approaching a constant. There is not really enough data to be certain about this empirical claim. (It can be established theoretically.) It seems clear that the convergence is neither linear nor quadratic (order = 1, 2), based on the rapid decrease of the fourth and fifth columns.
========================================================================= n p(n) e(n) e(n)/e(n-1) e(n)/e(n-1)^2 e(n)/e(n-1)^3 ========================================================================= 0 1.00000000 1.23606798 1 2.00000000 0.23606798 0.19098301 0.15450850 0.12500000 2 2.23529412 0.00077386 0.00327812 0.01388635 0.05882353 3 2.23606798 0.00000000 0.00000003 0.00003871 0.05002595 4 2.23606798 0.00000000 NaN NaN NaN =========================================================================