07-Jan-2022 19:16:38 ERRORS_TEST(): MATLAB/Octave version 9.8.0.1380330 (R2020a) Update 2 Test ERRORS. HOW RELIABLE ARE OUR CALCULATIONS? This program demonstrates how unreliable certain calculations are, when performed with a finite accuracy on a computer. Standard software is used where possible for the solution of these problems. Paraphrasing the reference: "Despite the spread of computers and calculators, many users suffer a certain computational gullibility. Whatever the computer produces is regarded as almost mathematically proven. It is often quite astonishing to discover that even in simple calculations with just a few operations it is possible to produce completely incorrect results, and that this is a necessary consequence of current computer techniques - not just the user's method, but the arithmetic processes of the computer as well. The following examples show how few operations are needed for a computer to return meaningless results. This fact applies not just to hand calculators but also to the biggest supercomputers. The worst aspect of these effects is that one never knows whether, when or where these sorts of errors will enter into a long, involved calculation and destroy the results. Only in the last few years have methods and solution processes been developed to mathematically describe these effects, and numerically control them." +=============================================+ | | | Exercise 1 | | | | P = 665857 | | Q = 470832 | | | | Compute: | | | | R = P^2 - 2 * Q^2 | | | | Correct value: | | | | R = 1 | | | +=============================================+ R = 1.0000000000000000 +=============================================+ | | | Exercise 2 | | | | P = 10864 | | Q = 18817 | | | | Compute: | | | | R = 9 * P^4 - Q^4 + 2 * Q^2 | | | | Correct value: | | | | R = 1 | | | +=============================================+ R = 2.0000000000000000 +=============================================+ | | | Exercise 3 | | | | P = 10 ^ 34 | | Q = - 2 | | | | Compute: | | | | R = ( P + Q ) - P | | | | Correct value: | | | | R = - 2 | | | +=============================================+ R = 0.0000000000000000 +=============================================+ | | | Exercise 4 | | | | Q = 1.091608 | | | | Compute: | | | | R = | | 170.4 * Q^3 | | - 356.41 * Q^2 | | + 168.97 * Q | | + 18.601 | | | | Correct value: | | | | R = 0.821248E-13 | | | +=============================================+ Direct evaluation R = 2.842170943040401e-14 Horner evaluation R = 6.394884621840902e-14 +=============================================+ | | | Exercise 5 | | | | Q = 0.707107 | | | | Compute: | | | | R = | | 8118 * Q^4 | | - 11482 * Q^3 | | + Q^2 | | + 5741 * Q | | - 2030 | | | | Correct value: | | | | R = - 0.191527325270E-10 | | | +=============================================+ Direct evaluation R = -1.932676241267473e-11 Horner evaluation R = -1.955413608811796e-11 +=============================================+ | | | Exercise 6 | | | | A = [ 64919121 - 159018721 ] | | [ 41869520.5 - 102558961 ] | | | | B = [ 1 ] | | [ 0 ] | | | | Solve: | | | | A * X = B | | | | Correct X: | | | | 205,117,922.0 | | 83,739,041.0 | | | +=============================================+ [Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.668288e-16.] [> In errors_test06 (line 51) In errors_test (line 34) In run (line 91) ] Exact x, computed x, residual: 1: 2.05118e+08 1.45867e+08 1 2: 8.3739e+07 5.955e+07 0 +=============================================+ | | | Exercise 6.5 | | | | A = [ 888445 887112 ] | | [ 887112 885781 ] | | | | det ( A ) = 1 | | | | B = [ 1 ] | | [ 0 ] | | | | Solve: | | | | A * X = B | | | | Correct X: | | | | 885781 | | -887112 | | | +=============================================+ Exact X, approx X, residual: 1: 885781 885805 -0.00012207 2: -887112 -887136 -0.00012207 +=============================================+ | | | Exercise 7 | | | | A = [ -367296 -43199 519436 -954302 ]| | [ 259718 -477151 -367295 -1043199 ]| | [ 886731 88897 -1254026 -1132096 ]| | [ 627013 566048 -886732 911103 ]| | | | B = [ 1 ] | | [ 1 ] | | [ 1 ] | | [ 0 ] | | | | Solve: | | | | A * X = B | | | | Correct X: | | | | 8.86731088897E+17 | | 8.86731088897E+11 | | 6.27013566048E+17 | | 6.27013566048E+11 | | | +=============================================+ [Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.556189e-18.] [> In errors_test07 (line 57) In errors_test (line 36) In run (line 91) ] x=A\b solution: 1: 1.90642e+10 2: 19064.2 3: 1.34804e+10 4: 13480.4 +=============================================+ | | | Exercise 8 | | | | P(X) = | | 170.4 * x^3 | | - 356.41 * x^2 | | + 168.97 * x | | + 18.601 | | | | Seek the minimizer X* of P(X) | | over the interval 0 <= X <= 2. | | | | The true minimizer lies in the interval | | | | [1.091607978, 1.091607981 ] | | | +=============================================+ FMINBND finds a minimizer at X = 1.091599046256172 with function value P(X) = 1.609067368235628e-08 n = 3 +=============================================+ | | | Exercise 9 | | | | P(X) = | | 223200658 * X^3 | | - 1083557822 * X^2 | | + 1753426039 * X | | - 945804881 | | | | Evaluate at 11 equally spaced values | | from 1.61801916 to 1.61801917 | | | | Correct table: | | | | ===== X ===== ========= P(X) ======== | | | | 1.618019160 - 0.17081105112320e-11 | | 1.618019161 - 0.89804011575510e-12 | | 1.618019162 - 0.34596536943057e-12 | | 1.618019163 - 0.51884933054474e-13 | | 1.618019164 - 0.15797467422848e-13 | | 1.618019165 - 0.23770163333175e-12 | | 1.618019166 - 0.71759609157723e-12 | | 1.618019167 - 0.14554795029553e-11 | | 1.618019168 - 0.24513505282621e-11 | | 1.618019169 - 0.37052078282936e-11 | | 1.618019170 - 0.52170500638460e-11 | | | +=============================================+ R1: compute P(X) the regular way. R2: compute P(X) using Horner's method. X, R1, R2 1.61801916 2.384185791015625e-07 -1.192092895507812e-07 1.618019161 0 -1.192092895507812e-07 1.618019162 0 1.192092895507812e-07 1.618019163 3.576278686523438e-07 -1.192092895507812e-07 1.618019164 -1.192092895507812e-07 1.192092895507812e-07 1.618019165 -1.192092895507812e-07 0 1.618019166 -1.192092895507812e-07 1.192092895507812e-07 1.618019167 -2.384185791015625e-07 0 1.618019168 -2.384185791015625e-07 1.192092895507812e-07 1.618019169 -3.576278686523438e-07 -1.192092895507812e-07 1.61801917 1.192092895507812e-07 -1.192092895507812e-07 +=============================================+ | | | Exercise 10 | | | | For the data: | | | | X Y | | | | 5201477 99999 | | 5201478 100000 | | 5201479 100001 | | | | Find the "nearest" straight line: | | | | Y(X) = A * X + B | | | | in the least squares sense. | | | | Correct values: | | | | A = 1.0 | | B = - 5101478 | | | | Evaluate: | | | | Y(5201480) | | | | Correct value: | | | | Y(5201480) = 100002.0 | | | +=============================================+ Slope A Intercept B Y(5201480) 1 -5101478 100002 +=============================================+ | | | Exercise 11 | | | | P(X) = | | 2124476931 * X^4 | | - 1226567328 * X^3 | | - 708158977 * X^2 | | + 408855776 * X | | + 1.0E-27 | | | | P(X) is positive for X >= 0. | | | | Use a minimizer to seek the minimum | | value of P(X) in (0.56, 0.59). | | | +=============================================+ FMINBND finds a minimizer at X = 0.5773583938327028 with function value P(X) = 0.09349286339359712 +=============================================+ | | | Exercise 12 | | | | Solve a linear system A * x = b | | | | A is the Hilbert matrix of order 9, | | multiplied by the least common multiple of | | 1 through 9, so all entries are integers. | | | | B = (1,0,0,...,0) | | | | We are interested in the first component | | X(1) of the solution vector. | | | | Correct X(1): | | | | 6.611036022800E-06 | | | +=============================================+ First solution component only: Exact X1(1) = 6.6110360228e-06 Computed X2(1) = 0.06428564667507462 +=============================================+ | | | Exercise 13 | | | | Solve a linear system A * x = b | | | | A is the Hilbert matrix of order 21 | | multiplied by the least common multiple of | | 1 through 25 so all entries are integral. | | | | B is (1,0,0,...,0) | | | | We are interested in the first component | | X(1) of the solution vector. | | | | Correct X(1): | | | | 2.013145339298E-15 | | | +=============================================+ [Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.605026e-19.] [> In errors_test13 (line 57) In errors_test (line 43) In run (line 91) ] First solution component only: Exact X1(1) = 2.013145339298e-15 Computed X2(1) = 3.040318275136052e-06 +=============================================+ | | | Exercise 14 | | | | A = [ 64079 57314 ] | | [ 51860 46385 ] | | | | B = [ 2 ] | | [ 305 ] | | | | Solve: | | | | A * X = B | | | | Correct X: | | | | - 46368.0 | | 51841.0 | | | +=============================================+ Exact X, computed X, Residual 1: -46368 -46368 0 2: 51841 51841 0 +=============================================+ | | | Exercise 15 | | | | F(X) = ( 4970 * X - 4923 ) | | -------------------------------- | | ( 4970 * X^2 - 9799 * X + 4830 ) | | | | Approximate F"(X) at X = 1 by: | | | | Del2(X)(F,H) = | | ( F(X-H) - 2*F(X) + F(X+H) ) / H^2 | | | | Correct values: | | | | F"(1.0) = 94.0 | | | | Del2(1.0)(F,1.0E-4) = 70.78819 | | Del2(1.0)(F,1.0E-5) = 93.76790 | | Del2(1.0)(F,1.0E-8) = 94.0 | | | +=============================================+ H Del2(X,H) 0.0001 70.78724237885581 1e-05 93.16032389961036 1e-08 894715.412869118 +=============================================+ | | | Exercise 16 | | | | P(X,Y) = | | 83521 * y^8 | | + 578 * x^2 * y^4 | | - 2 * x^4 | | + 2 * x^6 | | - x^8 | | | | Evaluate P(X,Y) for | | | | X = 9478657 | | Y = 2298912 | | | | Correct value: | | | | P(X,Y) = - 179,689,877,047,297 | | | +=============================================+ P(X,Y) = -1.088903574147003e+40 Use a scale factor on X and Y of 1e+06 Scaled value of P(X,Y) = 1.450476988585908e+18 Computed value of P(X,Y) = 1.450476988585908e+18 * 1e+06 ^ 8 +=============================================+ | | | Exercise 17 | | | | A = | | | | 2.718281828 | | -3.141592654 | | 1.414213562 | | 0.5772156649 | | 0.3010299957 | | | | B = | | | | 1486.2497 | | 878366.9879 | | - 22.37492 | | 4773714.647 | | 0.000185049 | | | | Compute the scalar product: | | | | R = A dot B | | | | Correct value: | | | | R = - 0.100657107E-10 | | | +=============================================+ R1 = standard loop = 1.025188136829667e-10 R2 = sum(a(1:n) .* b(1:n) ) = 1.025188136829667e-10 R3 = A' * B = 1.025188136829667e-10 +=============================================+ | | | Exercise 18 | | | | X = 192119201 | | Y = 35675640 | | | | Z(X,Y) = | | ( | | 1682 * x * y^4 | | + 3 * x^3 | | + 29 * x * y^2 | | - 2 * x^5 | | + 832 | | ) / 107751 | | | | Correct Z: | | | | 1783.0 | | | +=============================================+ Z = 0.00772150606490891 We use a scale factor on X and Y of 1000000 Scaled value of Z = 7.721506064908911e-33 Computed value of Z = 7.721506064908911e-33 * 1000000 ^ 5 +=============================================+ | | | Exercise 19 | | | | F1(X) = 1 - Cos(X) | | F2(X) = Sin^2(X) / ( 1 + Cos(X) ) | | F3(X) = Taylor series for 1 - Cos(X) | | around X = 0, up to X^6 terms. | | | | For all X | | | | F1(X) = F2(X) | | | | For X near 0, | | | | F3(X) should approximate F1(X) | | | +=============================================+ X 1-Cos Sin^2/(1+Cos) Taylor 0.5 0.1224174381096272 0.1224174381096273 0.1224175347222222 0.125 0.007802332770670994 0.007802332770670948 0.00780233277214898 0.03125 0.0004882415148635966 0.0004882415148636308 0.0004882415148636533 0.0078125 3.051742290494097e-05 3.05174229048867e-05 3.051742290488669e-05 0.001953125 1.907348026519706e-06 1.907348026482776e-06 1.907348026482776e-06 0.00048828125 1.192092872193129e-07 1.192092871823055e-07 1.192092871823055e-07 0.0001220703125 7.450580596923828e-09 7.45058058767197e-09 7.450580587671969e-09 3.0517578125e-05 4.656612873077393e-10 4.656612872715992e-10 4.656612872715992e-10 7.62939453125e-06 2.91038304567337e-11 2.910383045659253e-11 2.910383045659253e-11 1.9073486328125e-06 1.818989403545856e-12 1.818989403545305e-12 1.818989403545305e-12 4.76837158203125e-07 1.13686837721616e-13 1.136868377216139e-13 1.136868377216139e-13 1.192092895507812e-07 7.105427357601002e-15 7.105427357600994e-15 7.105427357600993e-15 2.980232238769531e-08 4.440892098500626e-16 4.440892098500626e-16 4.440892098500626e-16 +=============================================+ | | | Exercise 20 | | | | P(X) = | | 0.00005 * x^2 | | + 100 * x | | + 0.005 | | | | Evaluate the standard quadratic formula: | | | | X1 = 0.5 * ( -B + sqrt(B^2-4*A*C) ) / A | | X2 = 0.5 * ( -B - sqrt(B^2-4*A*C) ) / A | | | | Evaluate an alternate quadratic formula: | | | | X3 = 2 * C / ( - B + sqrt(B^2-4*A*C) ) | | X4 = 2 * C / ( - B - sqrt(B^2-4*A*C) ) | | | | Correct results: | | | | Root1 = -0.00005... | | Root2 = -1999999.99995... | | | +=============================================+ Standard quadratic formula: R1 = -5.000003966415534e-05, R2 = -1999999.99995 Alternate quadratic formula: R1 = -5.000000000125e-05, R2 = -1999999.99995 +=============================================+ | | | Exercise 21 | | | | P(X) = | | 1 * x^2 | | - 4 * x | | + 3.9999999 | | | | Evaluate the standard quadratic formula: | | | | X1 = 0.5 * ( -B + sqrt(B^2-4*A*C) ) / A | | X2 = 0.5 * ( -B - sqrt(B^2-4*A*C) ) / A | | | | Evaluate an alternate quadratic formula: | | | | X3 = 2 * C / ( - B + sqrt(B^2-4*A*C) ) | | X4 = 2 * C / ( - B - sqrt(B^2-4*A*C) ) | | | | Correct results: | | | | 1.999683772 | | 2.000316228 | | | +=============================================+ Standard quadratic formula: R1 = 2.000316227765758, R2 = 1.999683772234242 Alternate quadratic formula: R1 = 1.999683772234242, R2 = 2.000316227765758 +=============================================+ | | | Exercise 22 | | | | Define F(X,N) to be the difference between | | EXP(X) and its Taylor series, multiplied | | by N!. | | | | Formula 1: | | | | F(X,N) = N! * | | ( EXP(X) - (1+X+X^2/2+...+X^N/N!) ) | | | | But we can also determine a second | | definition, namely: | | | | Formula 2: | | | | F(X,N) = N*F(X,N-1) - X^N | | | | Compare formula 1 with formula 2 at X = 1 | | | +=============================================+ N Formula 1 Formula 2 x^n/n 0 1.718281828459046 1.718281828459046 0 1 0.7182818284590455 0.7182818284590455 1 2 0.4365636569180911 0.4365636569180911 0.5 3 0.3096909707542741 0.3096909707542732 0.3333333333333333 4 0.2387638830170999 0.2387638830170928 0.25 5 0.1938194150855033 0.1938194150854642 0.2 6 0.1629164905128633 0.1629164905127851 0.1666666666666667 7 0.1404154335889274 0.1404154335894958 0.1428571428571428 8 0.1233234687080653 0.1233234687159666 0.125 9 0.1099112184408568 0.109911218443699 0.1111111111111111 10 0.09911218431568614 0.09911218443698999 0.1 11 0.09023402737966535 0.09023402880688991 0.09090909090909091 12 0.08280827528324153 0.08280834568267892 0.08333333333333333 13 0.07650624909274484 0.07650849387482594 0.07692307692307693 14 0.07108062700353912 0.07111891424756323 0.07142857142857142 15 0.0662025447581982 0.06678371371344838 0.06666666666666667 16 0.0557495113753248 0.06853941941517405 0.0625 17 0 0.1651701300579589 0.05882352941176471 18 0 1.97306234104326 0.05555555555555555 19 0 36.48818447982194 0.05263157894736842 20 0 728.7636895964388 0.05 +=============================================+ | | | Exercise 23 | | | | Estimate EXP(-5.5) in two ways: | | | | R1(N) = sum ( I = 0 to N ) X^I/I | | R2(N) = 1 / sum ( I = 0 to N ) (-X)^I/I | | | | Correct value: | | | | EXP(-5.5) = 0.00408677143846406699 | | | +=============================================+ R1 = 0.004086771438461079 R2 = 0.004086771438464068 +=============================================+ | | | Exercise 24 | | | | Estimate the matrix exponential e^A. | | | | The matrix A is | | | | -49 24 | | -64 31 | | | | Correct value of e^A: | | | | -0.735759 0.551819 | | -1.471518 1.103638 | | | +=============================================+ EXPM1(A): Col: 1 2 Row 1 : -0.735759 0.551819 2 : -1.47152 1.10364 EXPM2(A): Col: 1 2 Row 1 : -0.735759 0.551819 2 : -1.47152 1.10364 EXPM3(A): Col: 1 2 Row 1 : -0.735759 0.551819 2 : -1.47152 1.10364 +=============================================+ | | | Exercise 26 | | | | Estimate an integer power of a | | real number. | | | | The computation is | | | | (1.0000001)^(2^27) | | | | Correct value is: | | | | 674,530.4707 | | | +=============================================+ X^(2^N) = 674530.476027064 X*X*...*X = 674530.4755217875 E^(logX*27^N)) = 674530.4760270639 E^(logX*E^(27*logN)) = 674530.47602707 +=============================================+ | | | Exercise 27 | | | | TAN(X) and ATAN(X) are inverse functions. | | Therefore, it should always be that: | | | | TAN(ATAN(X)) = X | | | +=============================================+ X TAN(ATAN(X)) Error 1 1 1.11022e-16 10 10 1.06581e-14 100 100 1.00897e-12 1000 1000 5.4456e-11 10000 10000 9.46056e-09 100000 100000 1.59882e-07 1e+06 1e+06 2.0701e-05 1e+07 1e+07 0.0102425 1e+08 1e+08 0.00457631 1e+09 1e+09 78.0719 1e+10 9.99999e+09 6950.63 1e+11 9.99994e+10 620594 1e+12 1.00007e+12 7.19169e+07 1e+13 1.00019e+13 1.86989e+09 1e+14 9.94704e+13 5.29576e+11 1e+15 1.05328e+15 5.32849e+13 1e+16 1.63312e+16 6.33124e+15 1e+17 1.63312e+16 8.36688e+16 1e+18 1.63312e+16 9.83669e+17 1e+19 1.63312e+16 9.98367e+18 1e+20 1.63312e+16 9.99837e+19 +=============================================+ | | | Exercise 28 | | | | f(x,y)= 333.75 y^6 | | + x^2 (11x^2y^2-y^6-121y^4-2) | | + 5.5y^8 + x/2/y | | | | Evaluate f(x,y) at: | | | | x = 77617, y = 33096. | | | | The correct answer: | | | | f(77617,33096) = -0.827... | | | +=============================================+ Computed result = -1.18059e+21 Correct result (to 3 digits) = -0.827 ERRORS_TEST: Normal end of execution. 07-Jan-2022 19:16:38