Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 15
Текст из файла (страница 15)
In this section, we discuss some of the manypossible ways of terminating iteration in a reasonable way and givetranslations of Algorithms 3.1 and 3.3, into FORTRAN.FORTRAN SUBROUTINE FOR THE BISECTIONALGORITHM 3.1SUBROUTINE BISECT ( F, A, B, XTOL, IFLAG )C****** I N P U T ******C F NAME OF FUNCTION WHOSE ZERO IS SOUGHT. NAME MUST APPEAR IN ANCE X T E R N A L STATEMENT IN THE CALLING PROGRAM.C A,B ENDPOINTS OF THE INTERVAL WHEREIN A ZERO IS SOUGHT.C XTOL DESIRED LENGTH OF OUTPUT INTERVAL.C****** O U T P U T ******C A,B ENDPOINTS OF INTERVAL KNOWN TO CONTAIN A ZERO OF F .82THE SOLUTION OF NONLINEAR EQUATIONSC IFLAG AN INTEGER,= -1, FAILURE SINCE F HAS SAME SIGN AT INPUT POINTS A AND BC= 0 , TERMINATION SINCE ABS(A-B)/2 .LE.
XTOLCC= l , TERMINATION SINCE ABS(A-B)/2 IS SO SMALL THAT ADDITIONTO (A+B)/2 MAKES NO DIFFERENCE .Cc****** M E T H O D ******C THE BISECTION ALGORITHM 3.1 IS USED, IN WHICH THE INTERVAL KNOWN TOC CONTAIN A ZERO IS REPEATEDLY HALVED .INTEGER IFLAGERROR,FA,FM,XMREAL A,B,F,XTOL,FA = F(A)IF (FA+F(B) .GT. 0.) THENIFLAG = -1PRINT 601,A,BFORMAT('F(X) IS OF SAME SIGN AT THE TWO ENDPOINTS',2E15.7)601RETURNEND IFCCCCERROR = ABS(B-A)DO WHILE ERROR .GT. XTOLERROR = ERROR/2.6IF (ERROR .LE. XTOL)RETURNXM = (A+B)/2.CHECK FOR UNREASONABLE ERROR REQUIREMENTIF (XM + ERROR .EQ. XM) THENIFLAG = 1RETURNEND IFFM = F(XM)CHOOSE NEW INTERVALIF (FA*FM .GT. 0.) THENA = XMFA = FMELSEB = XMEND IFGO TO 6ENDThe following program makes use of this subroutine to find the root ofEq.
(3.2), discussed in the preceding section.CMAIN PROGRAM FOR TRYING OUT BISECTION ROUTINEINTEGER IFLAGREAL A,B,ERROR,XIEXTERNAL FFA = 1.B = 2.CALL BISECT ( FF, A, B, 1.E-6, IFLAG )IF (IFLAG .LT. 0)STOPXI = (A+B)/2.ERROR = ABS(A-B)/2.PRINT 600, XI,ERROR600 FORMAT(' THE ZERO IS ',E15.7,' PLUS/MINUS ',E15.7)STOPENDREAL FUNCTION FF(X)REAL XFF = -1. - X*(1. - X*X)PRINT 600,X,FF600 FORMAT(' X, F(X) = ',2E15.7)RETURNENDWe now comment in some detail on the subroutine BISECT above.We have dropped the subscripts used in Algorithm 3.1. At any stage, the3.2FORTRAN PROGRAMS FOR SOME ITERATIVE METHODS83variables A and B contain the current lower and upper bound for the rootto be found, the initial values being supplied by the calling program.
Inparticular, the midpointis always the current best estimate for the root, its absolute difference fromthe root always being bounded byIteration is terminated oncewhere XTOL is a given absolute error bound. The calling program thenuses the current value of A and B to estimate the root. In addition to A, Band XTOL, the calling program is also expected to supply the FORTRANname of the function f(x) whose zero is to be located. Since the assumptionthat f(A) and f(B) are of opposite sign is essential to the algorithm, there isan initial test for this condition. If f(A) and f(B) are not of opposite sign,the routine immediately terminates.
The output variable IFLAG is used tosignal this unhappy event to the calling program.The subroutine never evaluates the given function more than once forthe same argument, but rather saves those values which might be needed insubsequent steps. This is a reasonable policy since the routine might wellbe used for functions whose evaluation is quite costly.
Finally, the routinehas some protection against an unreasonable error requirement: Suppose,for simplicity, that all calculations are carried out in four-decimal-digitfloating-point arithmetic and that the bounds A and B have already beenimproved to the point thatso thatThendepending on how rounding to four decimal places is done. In any event,so that, at the end of this step, neither A nor B has changed. If now thegiven error tolerance XTOL were less than 0.05, then the routine wouldnever terminate, since |B - A|/2 would never decrease below 0.05.
Toavoid such an infinite loop due to an unreasonable error requirement84THE SOLUTION OF NONLINEAR EQUATIONS(unreasonable since it requires the bounds A and B to be closer togetherthan is possible for two floating-point numbers of that precision to bewithout coinciding), the routine calculates the current value of ERROR asfollows. Initially,At the beginning of each step, ERROR is then halved, since that is thereduction in error per step of the bisection method. The routine terminates,once ERROR is so small that its floating-point addition to the currentvalue of XM does not change XM.Next we consider the modified regula falsi algorithm 3.3. In contrast tothe bisection method, the modified regula falsi is not guaranteed toproduce as small an interval containing the root as is possible with thefinite-precision arithmetic used (see Exercise 3.2-l).
Hence additionaltermination criteria must be used for this algorithm.FORTRAN PROGRAM USING THE MODIFIED REGULAFALSI ALGORITHM 33SUBROUTINE MRGFLS ( F, A, B, XTOL, FTOL, NTOL, W, IFLAG )C****** I N P U T ******C F NAME OF FUNCTION WHOSE ZERO IS SOUGHT. NAME MUST APPEAR IN ANE X T E R N A L STATEMENT IN THE CALLING PROGRAM .CC A,B ENDPOINTS OF INTERVAL WHEREIN ZERO IS SOUGHT.C XTOL DESIRED LENGTH OF OUTPUT INTERVALC FTOL DESIRED SIZE OF F(W)C NTOL NO MORE THAN NTOL ITERATION STEPS WILL BE CARRIED OUT.C****** O U T P U T ******C A,B ENDPOINTS OF INTERVAL CONTAINING THE ZERO .C W BEST ESTIMATE OF THE ZERO .C IFLAG AN INTEGER,C=-1, FAILURE, SINCE F HAS SAME SIGN AT INPUT POINTS A, B .= 0, TERMINATION BECAUSE ABS(A-B) .LE.
XTOL .C= 1, TERMINATION BECAUSE ABS(F(W)) .LE. FTOL .C= 2, TERMINATION BECAUSE NTOL ITERATION STEPS WERE CARRIED OUT .CC****** M E T H O D ******C THE MODIFIED REGULA FALSI ALGORITHM 3.3 IS USED. THIS MEANS THAT,C AT EACH STEP, LINEAR INTERPOLATION BETWEEN THE POINTS (A, FA) ANDC (B ,FB) IS USED, WITH FA*FB .LT. 0 ,TO GET A NEW POINT (W,F(W))C WHICH REPLACES ONE OF THESE IN SUCH A WAY THAT AGAIN FA*FB .LT. 0.C IN ADDITION, THE ORDINATE OF A POINT STAYING IN THE GAME FOR MOREC THAN ONE STEP IS CUT IN HALF AT EACH SUBSEQUENT STEP.INTEGER IFLAG,NTOL, NREAL A,B,F,FTOL,W,XTOL,FA,FB,FW,SIGNFA,PRVSFWFA = F(A)SIGNFA = SIGN(1., FA)FB = F(B)IF (SIGNFA*FB .GT.
0.) THENPRINT 601,A,B601FORMAT(' F(X) IS OF SAME SIGN AT THE TWO ENDPO INTS' ,2E15.7)IFLAG = -1RETURNEND IFCW - AFW = FADO 20 N=l,NTOL3.2FORTRAN PROGRAMS FOR SOME ITERATIVE METHODS85CHECK IF INTERVAL IS SMALL ENOUGH.IF (ABS(A-B) .LE. XTOL) THENIFLAG = 0RETURNEND IFCHECK IF FUNCTION VALUE AT W IS SMALL ENOUGH .CIF (ABS(FW) .LE. FTOL) THENIFLAG = 1RETURNCEND IFGET NEW GUESS W BY LINEAR INTERPOLATION .W = (FA*B - FB*A)/(FA - FB)PRVSFW = SIGN(1.,FW)FW = F(W)CCHANGE TO NEW INTERVALIF (SIGNFA*FW .GT. 0.) THENA = WFA = FWIF (FW*PRVSFW .GT. 0.) FB = FB/2.ELSEB = WFB = FWIF (FW*PRVSFW .GT. 0.) FA = FA/2.END IFCONTINUEPRINT 620,NTOL620 FORMAT(' NO CONVERGENCE IN ',I5,' ITERATIONS')IFLAG = 2RETURNENDCFirst, the routine terminates if the newly computed function value isno bigger in absolute value than a given tolerance FTOL.
This brings inthe point of view that an “approximate root” of the equation f(x) = 0 is apoint x at which |f(x)| is “small.” Also, since the routine repeatedly dividesby function values, such a termination is necessary in order to avoid, inextreme cases, division by zero.Second, the routine terminates when more than a given number NTOLof iteration steps have been carried out. In a way, NTOL specifies theamount of computing users are willing to invest in solving their problems.Use of such a termination criterion also protects users against unreasonable error requirements and programming errors, and against the possibility that they have not fully understood the problem they are trying tosolve.
Hence such a termination criterion should be used with any iterativemethod.As in the routine for the bisection method, the subroutine MRGFLSreturns an integer IFLAG which indicates why iteration was terminated,and the latest value of the bounds A and B for the desired root. Finally, aswith the bisection routine, the routine never evaluates the given functionmore than once for the same argument.Algorithms 3.4 and 3.5 for the secant method and Newton’s method,respectively, do not necessarily bracket a root. Rather, both generate asequence x0, x1, x2, . . . , which, so one hopes, converges to the desiredroot of the given equation f(x) = 0. Hence both algorithms should beviewed primarily as finding points at which f(x) is “small” in absolute86THE SOLUTION OF NONLINEAR EQUATIONSvalue; iteration is terminated when the newly computed function value isabsolutely less than a given FTOL.The iteration may also be terminated when successive iterates differin absolute value by less than a given number XTOL.
It is customarytherefore to use one or both of the following termination criteria for eitherthe secant or Newton’s method:(3.13)If the size of the numbers involved is not known in advance, it is usuallybetter to use relative error requirements, i.e., to terminate if(3.14)where FSIZE is an estimate of the magnitude of f(x) in some vicinity ofthe root established during the iteration.In Sec. 1.6 we discussed the danger of concluding that a givensequence has “converged” just because two successive terms in thesequence differ by “very little.” Such a criterion is nevertheless commonlyused in routines for the secant and Newton methods. For one thing, such acriterion is necessary in the secant method to avoid division by zero.
Also,in both methods, the difference between the last two iterates calculated is arather conservative bound for the error in the most recent iterate once theiterates are “close enough” to the root. To put it naively: If successiveiterates do not differ by much, there is little reason to go on iterating.Subroutines for the Newton and secant methods are not included in thetext but are left as exercises for the student.Example 3.2a Find the real positive root of the equationThe results for Algorithms 3.1, 3.3, 3.4, and 3.5 are given in the following table, whichparallels the table in Example 3.1.3.2FORTRAN PROGRAMS FOR SOME ITERATIVE METHODS87Example 3.2b The so-called biasing problem in electronic circuit design requires thesolution of an equation of the formwhere v representsthe voltage, I is a measure of current, and q is a parameter relating the electron chargeand the absolute temperature. In a typical engineering problem this equation wouldneed to be solved for various values of the parameters I and q to see how the smallestpositive zero of f(v) changes as the parameters change.Using Newton’s method find the smallest positive zero of f(v) under two differentsets of parameter values (I, q) = (10-8, 40) and (I, q) = (10-6, 20).
Set XTOL = 10-8and FTOL = 10 -7 .The results using the indicated starting values are given below.In this example a poor selection of starting values will lead to divergence.EXERCISES3.2-l Try to find the root x = 1.3333 of the equation (x - 1.3333)3 = 0 to five places ofaccuracy using the modified regula falsi algorithm 3.3 and starting with the interval [1,2].Why does the method fail in this case to give a “small” interval containing the root?3.2-2 Because of the use of the product FA*FM in the subroutine BISECT, overflow orunderflow may occur during the execution of this subroutine, even though the function valuesFA and FM are well-defined floating-point numbers.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
