Thompson - Computing for Scientists and Engineers (523188), страница 23
Текст из файла (страница 23)
For example, if x1 = 1000000 and x2 = 999999, then their difference x1- x 2 = 1, has only one significant figure, whereas x1 and x2 had about sixsignificant figures. We now examine two examples in which such perils of subtractive-cancellation errors may arise, variance calculations and the solution of quadraticequations.4.3NUMERICAL NOISE IN COMPUTING117Variance calculations (and thence standard deviations) among large numbers provide good examples of difficulties from subtractive cancellation.
Suppose that theaverage of a set of x values is calculated as(4.14)The variance among this set, V, is defined by(4.15)Analytically, this can be expanded to produce(4.16)Exercise 4.7Verify the algebraic equivalence of these two expressions for the variance, V. nNotice that in (4.15) all N of the xj values must be input and their average must becalculated before their variance can be calculated, whereas in (4.16) the sum ofsquares and the sum of values (needed to compute the average) may be accumulatedas soon as each x value is entered. If (4.15) is used for large N, more storage thansuch devices as pocket calculators have available for data will often be required. So,for calculators the second method is almost always used, because it requires storageonly for the running sums.
Numerically, the outcome of using (4.16) can be complete loss of significance, as you will discover if you work the following exercise.Exercise 4.8(a) Use a pocket calculator to calculate the average and variance for the threedata values x1 = 999999, x2 = 1000000, x3 = 1000001. Show that using(4.16) you should get V = 2. What result does your calculator get, and why?(b) Reduce the order of magnitude of the three data by successive powers of 10until the two formulas agree.
From this numerical experiment estimate the number of significant figures carried by your calculator. Does this agree with its instruction manual? nBy using my pocket calculator and the data in Exercise 4.8 (a), I obtained a variance of zero according to (4.16). With values of 99999, 100000, and 100001,however, I obtained the correct variance by using either (4.15) or (4.16).
My calculator manual advertises 10 digits for entry or display and 12 digits for each step of acalculation: I now believe it, but I don’t like the subtractive cancellation effects. Anextensive discussion of roundoff errors in computers is given in Chapter 1 of thetextbook by Maron.118NUMERICAL DERIVATIVES AND INTEGRALSQuadratic equations provide another example where subtractive cancellation canbe a severe problem with a numerical calculation as simple as finding their roots accurately.
Write the quadratic equation with two real roots as(4.17)The choice of unity for the first root, x = x 1 = 1, can always be achieved by appropriate scaling of x, and is convenient when comparing the second root,x = x2 = with the first, especially when is small in magnitude. By expanding(4.17) into standard quadratic-equation form, we have(4.18)The usual way that this is solved is to write the solution as(4.19)Since we know the roots, we could simplify this analytically and obtain the two solutions given above. Proceeding numerically, however, is quite another story.
I used my pocket calculator to input then calculate the two numerical solutions in(4.19), xn1 and xn2, as written. This calculator enters and displays 10 decimal digitsbut it uses 12 digits at each step of a calculation. For values of between 10-l1 and10-5 (and presumably for larger values) the first root, xn1, was exact, but the relative error in the second root, xn2, behaved quite erratically, as Table 4.2 shows.TABLE 4.2 Relative errors in root of (4.18) calculated using (4.19) and a 10-digit calculator.This example illustrates that when is less than the square root of the accuracy ofthe calculator< 10-5 for a 10-digit calculator), then the smaller root cannot bereliably computed, because the squared factor under the root in (4.19) is notcalculated accurately enough. Significant digits are then lost when the root is subtracted from the first term in the denominator of (4.19).When solving quadratic equations, how can we reduce this numerical noise fromsubtractive cancellation? The easiest way is to avoid subtracting.
To show how thiscan be accomplished, and to cast the quadratic-equation problem in a more conventional form, let us write(4.20)4.3NUMERICAL NOISE IN COMPUTINGThe general analytical solution of this quadratic when a1190 is(4.2 1)Exercise 4.9(a) Multiply numerator and denominator of (4.21) by factororder to show that the roots of a quadratic are equivalently given byin(4.22)in which the denominator is assumed to be nonzero.(b) In order to make the derivation of (4.22) less ad hoc, consider the followingtransformation of the quadratic equation (4.20). Assume that the roots are nonzero (c 0), so that (4.20) can be divided throughout by x 2.
You now have aquadratic equation in the variable l/x with the roles of a and c interchanged.Write down the solution for this quadratic by using (4.21), then take its reciprocal. Viola! the solution (4.22) is obtained. nBetween the use of the two equations (4.21) and (4.22) for the roots, we can avoidsubtractive cancellation completely. For, if b > 0 (4.22) may be used with the +sign and (4.21) may be used with the - sign. If b < 0 the opposite choice of signsis appropriate.Program for roots of quadratic equationsThe program that we now provide uses the algorithm devised in Exercise 4.9. OurProgram 4.4, Quadratic Equation Roots, computes the roots of quadraticequations (4.20) for which the coefficients a, b, and c are each real numbers.Program Quadratic Equation Roots has a small driver program to controlprogram execution, and function quadroots to solve for the roots.
The main program reads a, b, and c; then if any one of them is nonzero it attempts a solution.Most of the function quadroots is occupied with taking care of the case a = 0,including the subcase that gives rise to inconsistent solutions, namely b = 0. Inthis subcase if the coefficients of the powers of x in (4.20) are both zero, then theroots are indeterminate if c = 0 and the equation is inconsistent if c 0.If a is nonzero, then subtractive cancellation can be a problem only if the squareroot is real, since there is no subtraction when the real and imaginary parts of a complex number have opposite signs, as discussed in Section 2.1.
Therefore, a test fora positive discriminant is made before deciding whether to use formulas (4.21) and(4.22) for the case of real roots. For complex roots the usual form of the solution,(4.2 l), is used for the two roots, which are related through complex conjugation.120NUMERICAL DERIVATIVES AND INTEGRALSPROGRAM 4.4 Quadratic equation roots avoiding subtractive cancellation.#include#include<stdio.h><math.h>main()/* Quadratic Equation Roots */doublea,b,c,ReXp,ImXp,ReXm,ImXm;int error;void quadroots();printf("Quadratic Equation Roots\n");a = 1;while ( a != 0 || b != 0 || c!= 0 )printf("\nInput a,b,c (all zero to end) :\n");scanf("%lf%lf%lf",&a,&b,&c);if ( a == 0 && b == 0 && c == 0 ){printf("\nEnd Quadratic Equation Roots"); exit(O);}quadroots(a,b,c,%ReXp,&ImXp,&ReXm,&ImXm,&error);printf("\nRoots are:");printf("(%lf)+i(%lf) & (%lf)+i(%lf)",ReXp,ImXp,ReXm,ImXm);printf("\nError flag = %i",error);void quadroots(a,b,c,ReXp,ImXp,ReXm,ImXm,error)/* Quadratic equation roots; Real and Imginary */doublea,b,c,*ReXp,*ImXp,*ReXm,*ImXm;int *error;{double Disc;*ReXp= *ImXp = *ReXm = *ImXm =0;*error = 0; /* no problems SO Far */if ( a == 0 ){if (b == 0 )if ( c != 0 ) *error = 2; /* inconsistent equation */*error = 1; /* roots indeterminate */else4.3NUMERICAL NOISE IN COMPUTING121else /* a is zero, b is not zero */{*ReXp = -c/b; *ReXm = *ReXp; /* degenerate */}}1else /* a is not zero */Disc = b*b-4*a*c; /* discriminant */if ( Disc >= 0 ) /* real roots */{if (b>= 0){*ReXp = -2*c/(b+sqrt(Disc));*ReXm = (-b-sqrt(Disc))/(2*a);else /* b is negative */*ReXp = (-bsqrt (Disc))/(2*a);*ReXm = -2*c/ (b-sqrt(Disc));else /* complex roots */&ReXp = -b/(2*a);*ImXp = sqrt(-Disc)/(2*a);*ReXm = *ReXp;*ImXm = *ImXp;After function quadroots has been executed, the two complex numbersReXp + i ImXp and ReXm + i ImXm are returned, together with the integer variable error, which is 0 for no problem, 1 for indeterminate roots, and 2 for an inconsistent equation.
This value is output with the roots before the program requestsanother set of coefficients. If all three coefficients are input as zero, then executionof Quadratic Equation Roots terminates.Exercise 4.10(a) Write down the solution of (4.20) when a = 0, and also when both a and bare zero. Check that these solutions agree with what is coded in the functionquadroots in the program Quadratic Equation Roots.(b) Test the program with some simple (say integer) values of a, b, and c, including the special cases of a = 0, as well as when both a and b are zero.122NUMERICAL DERIVATIVES AND INTEGRALS(c) Verify that subtractive cancellation is not a problem with the present algorithm. To do this, choose in Quadratic Equation Roots values of a = 1,b = -(1 + c =where << 1, as in the examples in Table 4.2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
