Chapter 3 Here is a VBA Implementation of the Algorithm

CHAPTER 3 3.1 Here is a VBA implementation of the algorithm: Option Explicit Sub GetEps() Dim epsilon As Single epsilon = 1 Do If epsilon + 1 <= 1 Then Exit Do epsilon = epsilon / 2 Loop epsilon = 2 * epsilon MsgBox epsilon End Sub It yields a result of 1.19209×10−7 on my desktop PC. 3.2 Here is a VBA implementation of the algorithm: Option Explicit Sub GetMin() Dim x As Single, xmin As Single x = 1 Do If x <= 0 Then Exit Do xmin = x x = x / 2 Loop MsgBox xmin End Sub It yields a result of 1.4013×10−45 on my desktop PC. 3.3 The maximum negative value of the exponent for a computer that uses e bits to store the exponent is emin = −(2e −1 − 1) Because of normalization, the minimum mantissa is 1/b = 2−1 = 0.5. Therefore, the minimum number is xmin = 2 −12 − ( 2 e−1 −1) = 2− 2 e−1 For example, for an 8-bit exponent xmin = 2 − 2 8−1 = 2 −128 = 2.939 × 10 − 39 This result contradicts the value from Prob. 3.2 (1.4013×10−45). This amounts to an additional 21 divisions (i.e., 21 orders of magnitude lower in base 2). I do not know the reason for the discrepancy. However, the problem illustrates the value of determining such quantities via a program rather than relying on theoretical values. 3.4 VBA Program to compute in ascending order Option Explicit Sub Series() Dim i As Integer, n As Integer Dim sum As Single, pi As Single pi = 4 * Atn(1) sum = 0 n = 10000 For i = 1 To n sum = sum + 1 / i ^ 2 Next i MsgBox sum 'Display true percent relatve error MsgBox Abs(sum - pi ^ 2 / 6) / (pi ^ 2 / 6) End Sub This yields a result of 1.644725 with a true relative error of 6.086×10−5. VBA Program to compute in descending order: Option Explicit Sub Series() Dim i As Integer, n As Integer Dim sum As Single, pi As Single pi = 4 * Atn(1) sum = 0 n = 10000 For i = n To 1 Step -1 sum = sum + 1 / i ^ 2 Next i MsgBox sum 'Display true percent relatve error MsgBox Abs(sum - pi ^ 2 / 6) / (pi ^ 2 / 6) End Sub This yields a result of 1.644725 with a true relative error of 1.270×10−4 The latter version yields a superior result because summing in descending order mitigates the roundoff error that occurs when adding a large and small number. 3.5 Remember that the machine epsilon is related to the number of significant digits by Eq. 3.11 ξ = b1−t which can be solved for base 10 and for a machine epsilon of 1.19209×10−7 for t = 1 − log10 (ξ ) = 1 − log10 (1.19209 ×10-7 ) = 7.92 To be conservative, assume that 7 significant figures is good enough. Recall that Eq. 3.7 can then be used to estimate a stopping criterion, ε s = (0.5 × 102 − n )% Thus, for 7 significant digits, the result would be ε s = (0.5 × 102 − 7 )% = 5 ×10 −6% The total calculation can be expressed in one formula as ε s = (0.5 × 102 − Int (1− log10 (ξ )) )% It should be noted that iterating to the machine precision is often overkill. Consequently, many applications use the old engineering rule of thumb that you should iterate to 3 significant digits or better. As an application, I used Excel to evaluate the second series from Prob. 3.6. The results are: Notice how after summing 27 terms, the result is correct to 7 significant figures. At this point, both the true and the approximate percent relative errors are at 6.16×10−6 %. At this point, the process would repeat one more time so that the error estimates would fall below the precalculated stopping criterion of 5×10−6 %. 3.6 For the first series, after 25 terms are summed, the result is The results are oscillating. If carried out further to n = 39, the series will eventually converge to within 7 significant digits. In contrast the second series converges faster. It attains 7 significant digits at n = 28. 3.9 Solution: 21 x 21 x 120 = 52920 words @ 64 bits/word = 8 bytes/word 52920 words @ 8 bytes/word = 423360 bytes 423360 bytes / 1024 bytes/kilobyte = 413.4 kilobytes = 0.41 M bytes 3.10 Solution: % % % Given: Taylor Series Approximation for cos(x) = 1 - x^2/2! + x^4/4! - ... Find: number of terms needed to represent cos(x) to 8 significant figures at the point where: x=0.2 pi x=0.2*pi; es=0.5e-08; %approximation cos=1; j=1; % j=terms counter fprintf('j= %2.0f cos(x)= %0.10f\n', j,cos) fact=1; for i=2:2:100 j=j+1; fact=fact*i*(i-1); cosn=cos+((-1)^(j+1))*((x)^i)/fact; ea=abs((cosn-cos)/cosn); if ea< es fprintf('j= %2.0f cos(x)= %0.10f ea = %0.1e CONVERGENCE es= %0.1e',j,cosn,ea,es) break end fprintf( 'j= %2.0f cos(x)= %0.10f ea = %0.1e\n',j,cosn,ea ) cos=cosn; end j= j= j= j= j= j= 1 2 3 4 5 6 cos(x)= cos(x)= cos(x)= cos(x)= cos(x)= cos(x)= 1.0000000000 0.8026079120 0.8091018514 0.8090163946 0.8090169970 0.8090169944 ea ea ea ea ea = = = = = 2.5e-001 8.0e-003 1.1e-004 7.4e-007 3.3e-009 CONVERGENCE es = 5.0e-009»