Applied Numerical Methods Using MATLAB. Won Y. Yang

      (1.2 17)

      (1.2 18)

      This implies that, in the worst‐case, the relative error in multiplication/division may be as large as the sum of the relative errors of the two numbers.

1.2.5 Tips for Avoiding Large Errors

In this section, we will look over several tips to reduce the chance of large errors occurring in calculations.

      First, to decrease the magnitude of round‐off errors and to lower the possibility of overflow/underflow errors, make the intermediate result as close to 1 (one) as possible in consecutive multiplication/division processes. According to this rule, when computing xy/z, we program the formula as

       %nm125_1.m x=36; y=1e16; for n=[-20 -19 19 20] fprintf('ŷ%2d/ê%2dx=%25.15e\n',n,n,ŷn/exp(n*x)); fprintf('(y/êx)̂%2d=%25.15e\n',n,(y/exp(x))̂n); end

      For instance, when computing yⁿ/e^nx with x > 1 and y > 1, we would program it as (y/e^x)ⁿ rather than as yⁿ/e^nx, so that overflow/underflow can be avoided. You may verify this by running the above MATLAB script “nm125_1.m”.

      >nm125_1 ŷ-20/ê-20x= 0.000000000000000e+000 (y/êx)̂-20= 4.920700930263814e-008 ŷ-19/ê-19x= 1.141367814854768e-007 (y/êx)̂-19= 1.141367814854769e-007 ŷ19/ê19x= 8.761417546430845e+006 (y/êx)̂19= 8.761417546430843e+006 ŷ20/ê20x= NaN (y/êx)̂20= 2.032230802424294e+007

      Second, in order to prevent ‘loss of significance’, it is important to avoid a ‘bad subtraction’ (Section 1.2.2), i.e. a subtraction of a number from another number having almost equal value. Let us consider a simple problem of finding the roots of a second‐degree equation ax² + bx + c = 0 by using the quadratic formula

(1.2 19)

      Let |4ac| ≪ b². Then, depending on the sign of b, a ‘bad subtraction’ may be encountered when we try to find x₁ or x₂, which is a smaller one of the two roots. This implies that it is safe from the ‘loss of significance’ to compute the root having the larger absolute value first and then obtain the other root by using the relation (between the roots and the coefficients) x₁x₂ = c/a (see Problem 1.20(b)).

       %nm125_2.m: roundoff error test f1=@(x)(1-cos(x))/x/x; % Eq.(1.2.20-1) f2=@(x)sin(x)*sin(x)/x/x/(1+cos(x)); % Eq.(1.2.20-2) for k=0:1 x=k*pi; tmp= 1; for k1=1:8 tmp=tmp*0.1; x1= x+tmp; fprintf('At x=%10.8f, ', x1) fprintf('f1(x)=%18.12e; f2(x)=%18.12e', f1(x1),f2(x1)); end end

For another instance, we consider the following two formulas, which are analytically the same, but numerically different.

      (1.2 20)

      It is safe to use f₁(x) for x ≈ π since the term (1 + cos x) in f₂(x) is a ‘bad subtraction’, while it is safe to use f₂(x) for x ≈ 0 since the term (1 − cos x) in f₁(x) is a ‘bad subtraction’. Let us run the above MATLAB script “nm125_2.m” to confirm this. Under is the running result. This implies that we might use some formulas to avoid a ‘bad subtraction’.

      >nm125_2 At x=0.10000000, f1(x)=4.995834721974e-01; f2(x)=4.995834721974e-01 At x=0.01000000, f1(x)=4.999958333474e-01; f2(x)=4.999958333472e-01 At x=0.00100000, f1(x)=4.999999583255e-01; f2(x)=4.999999583333e-01 At x=0.00010000, f1(x)=4.999999969613e-01; f2(x)=4.999999995833e-01 At x=0.00001000, f1(x)=5.000000413702e-01; f2(x)=4.999999999958e-01 At x=0.00000100, f1(x)=5.000444502912e-01; f2(x)=5.000000000000e-01 At x=0.00000010, f1(x)=4.996003610813e-01; f2(x)=5.000000000000e-01 At x=0.00000001, f1(x)=0.000000000000e+00; f2(x)=5.000000000000e-01 At x=3.24159265, f1(x)=1.898571371550e-01; f2(x)=1.898571371550e-01 At x=3.15159265, f1(x)=2.013534055392e-01; f2(x)=2.013534055391e-01 At x=3.14259265, f1(x)=2.025133720884e-01; f2(x)=2.025133720914e-01 At x=3.14169265, f1(x)=2.026294667803e-01; f2(x)=2.026294678432e-01 At x=3.14160265, f1(x)=2.026410772244e-01; f2(x)=2.026410604538e-01 At x=3.14159365, f1(x)=2.026422382785e-01; f2(x)=2.026242248740e-01 At x=3.14159275, f1(x)=2.026423543841e-01; f2(x)=2.028044503269e-01 At x=3.14159266, f1(x)=2.026423659946e-01; f2(x)= Inf

      It may be helpful for avoiding a ‘bad subtraction’ to use the Taylor series expansion [W-5] rather than using the exponential function directly for the computation of e^x. For example, suppose we want to find

(1.2 21)

      We can use the Taylor series expansion up to just the fourth‐order of e^x about x = 0: