Читать книгу Applied Numerical Methods Using MATLAB - Won Y. Yang - Страница 33

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

 (xy)/z when x and y in the multiplication are very different in magnitude,

 x(y/z) when y and z in the division are close in magnitude, and

 (x/z)y when x and z in the division are close in magnitude.

%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:

to approximate the above function (1.2.19) as

(1.2 22)

Noting that the true value of (1.2.21) is computed to be 1 by using the L'Hopital's rule [W-7], we run the MATLAB script “nm125_3.m” to find which one of the two formulas f₃(x) and f₄(x) is better for finding the value of the 1.2.21 at x = 0. Would you compare them based on the running result shown below? How can the approximate formula f₄(x) outrun the true one f₃(x) for the numerical purpose, though not usual? It is because the zero factors in the numerator/denominator of f₃(x) are cancelled to set f₄(x) free from the terror of a ‘bad subtraction’.

>nm125_3 At x=0.100000000000, f3(x)=1.051709180756e+00; f4(x)=1.084166666667e+00 At x=0.010000000000, f3(x)=1.005016708417e+00; f4(x)=1.008341666667e+00 At x=0.001000000000, f3(x)=1.000500166708e+00; f4(x)=1.000833416667e+00 At x=0.000100000000, f3(x)=1.000050001667e+00; f4(x)=1.000083334167e+00 At x=0.000010000000, f3(x)=1.000005000007e+00; f4(x)=1.000008333342e+00 At x=0.000001000000, f3(x)=1.000000499962e+00; f4(x)=1.000000833333e+00 At x=0.000000100000, f3(x)=1.000000049434e+00; f4(x)=1.000000083333e+00 At x=0.000000010000, f3(x)=9.999999939225e-01; f4(x)=1.000000008333e+00 At x=0.000000001000, f3(x)=1.000000082740e+00; f4(x)=1.000000000833e+00 At x=0.000000000100, f3(x)=1.000000082740e+00; f4(x)=1.000000000083e+00 At x=0.000000000010, f3(x)=1.000000082740e+00; f4(x)=1.000000000008e+00 At x=0.000000000001, f3(x)=1.000088900582e+00; f4(x)=1.000000000001e+00

%nm125_3.m: reduce the roundoff error using Taylor series f3=@(x)(exp(x)-1)/x; % LHS of Eq.(1.2.22) f4=@(x)((x/4+1)*x/3)+x/2+1; % RHS of Eq.(1.2.22) x=0; tmp=1; for k1=1:12 tmp=tmp*0.1; x1=x+tmp; fprintf('At x=%14.12f, ', x1) fprintf('f3(x)=%18.12e; f4(x)=%18.12e', f3(x1),f4(x1)); end