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 ex. For example, suppose we want to find
(1.2 21)
We can use the Taylor series expansion up to just the fourth‐order of ex 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 f3(x) and f4(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 f4(x) outrun the true one f3(x) for the numerical purpose, though not usual? It is because the zero factors in the numerator/denominator of f3(x) are cancelled to set f4(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
1.3 Toward Good Program
Among the various criteria about the quality of a general program, the most important one is how robust its performance is against the change of the problem properties and the initial values. A good program guides the program users who do not know much about the program and at least give them a warning message without runtime error for their minor mistake. There are many other features that need to be considered, such as user friendliness, compactness and elegance, readability. But, as far as the numerical methods are concerned, the accuracy of solution, execution speed (time efficiency), and memory utilization (space efficiency) are of utmost concern. Since some tips to achieve the accuracy or at least to avoid large errors (including overflow/underflow) are given in the last section, we will look over the issues of execution speed and memory utilization.
1.3.1 Nested Computing for Computational Efficiency
The execution speed of a program for a numerical solution depends mostly on the number of function (subroutine) calls and arithmetic operations performed in the program. Therefore, we like the algorithm requiring fewer function calls and arithmetic operations. For instance, suppose we want to evaluate the value of a polynomial
(1.3. 1)
It is better to use the nested (computing) structure (as follows) than to use the above form as it is
(1.3.2)
Note that the numbers of multiplications needed in Eqs. (1.3.2) and (1.3.1) are 4 and (4 + 3 + 2 + 1 = 9), respectively. This point is illustrated by the script “nm131_1.m”, where a polynomial polyval()’. Interested readers could run this script to see that the nested multiplication (Eq. (1.3.2)) and ‘ polyval()’ (fabricated in a nested structure) take less time than the plain method (Eq. (1.1)), while ‘ polyval()’ may take longer time because of some overhead time for being called.
%nm131_1.m: nested multiplication vs. plain multiple multiplication N=1000000+1; a=[1:N]; x=1; tic % initialize the timer p=sum(a.*x.̂[N-1:-1:0]); % Plain multiplication p time_plain=toc % Operation time for the sum of multiplications tic, pn=a(1); for i=2:N % Nested multiplication pn = pn*x +a(i); end pn, time_nested=toc % Operation time for the nested multiplications tic, polyval(a,x) time_polyval=toc % Operation time for using polyval()
Programming in a nested structure is not only recommended for time‐efficient computation, but also may be critical to the solution. For instance, consider a problem of finding the value
(1.3.3)
%nm131_2_1.m: nested structure lam=100; K=155; p=exp(-lam); tic S=0; for k=1:K p=p*lam/k; S=S+p; end S time_nested=toc
|
%nm131_2_2.m: not nested structure lam=100;
|
