Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 31
Текст из файла (страница 31)
Ifmore appear to be necessary, FLAG is set to 1 and the code terminates before theconvergence test is satisfied. The value FLAG = -1 indicates invalid input, that is,ABSERR < 0 or RELERR < 10u, and FLAG = -2 means F(B)F(C) > 0. The valueRESIDL (or residual in C and C++) is the final residual F(B). Convergence is judgedby the length of an interval known to contain a root. The algorithm is so robust thatit can locate roots of functions that are only piecewise continuous.
If it is applied toa function that has a pole of odd multiplicity, it might locate a pole rather than a root.This is recognized by a “large” residual and signaled by FLAG = 2.Example 4.4. The function F(x) = e-x - 2x has F(0) > 0 and F(1) < 0, hence theequation F(x) = 0 has a root between C = 0 and B = 1.
The sample program providedillustrates the use of the zero-finding routine. Note that a globally defined integer isused to count the number of F evaluations required. In FORTRAN this is accomplished via a COMMON statement (preferably labeled COMMON) and in C or C++it is done by an appropriate placement of the variable declarations. The output isExample 4.5.In Chapter 5 we discuss a problem that requires the solution offor its smallest positive root. The particular values of the parameters used there arex = y = z = 50 and a = l, b = 2, c= 100. Figure 4.5 is a rough sketch of φ for a2 <154CHAPTER 4ROOTS OF NONLINEAR EQUATIONSb2 < c2.
It is drawn by looking at the behavior of φ as λandThe portion of interest to us is between -a2 andSince asfrom theright, φ(λ)and since asthe continuous function φ musthave a zero larger than -a2. Differentiating φ(λ) we getBecause φ(λ) is strictly decreasing, there is only one zero λ0 greater than -a2. Theroot is simple sinceInterestingly, φ(λ) can be scaled to a cubic polynomialThis allows us to apply some bounds on the relative error of an approximate zero of apolynomial developed in Exercise 4.31.The equation φ(λ) = 0 was solved using the code Zero with relative and absoluteerror requests of 10-6 and 10-8, respectively.
Poor initial values of B and C were usedto show the excellent rate of convergence. Table 4.1 displays successive values of Band C and tells which method was used by Zero in computing B.Table 4.1. Solution of φ(λ) = 0 by ZeroFor the error bounds in Exercise 4.31, we computed P(B) = 1.6261594 × 104 andP´(B) = 8.5162684 × 107. The constant term in P(λ) is a0 = x2b2c2 - y2a2c2 z2a2b2 + a2b2c2, which in this case is - 1.2497000 × 108. The bounds state that thereis a root rj of P such thatand a root ri such thatThe second error bound is quite good, but the first is pessimistic. Generally we donot know which of the two bounds will be better. An interesting point here is the sizeof the residual P(B). The value B is supposed to be a root of φ(λ) = 0 and P(λ) =0.
If its quality were judged by the size of the residual, B might be thought a poor4.3 ROUTINES FOR ZERO FINDING155Figure 4.5 φ(λ) for Example 4.5.approximation when solving the one equation and a good approximation when solvingthe other. This is despite the fact that it approximates the same root in both cases. Tobe more concrete, MATHCAD’S default convergence test that the residual be smallerin magnitude than 10-3 would reject an approximate root of P(λ) = 0 that is actuallyquite accurate. This problem illustrates well the issue of scaling. The convergence testin Zero and in similar codes is reliable. The magnitude of the residual reported by thecode is important only for detecting the convergence to a pole that is possible whenthe function input is not continuous.
Otherwise its size is a statement about the scalingof the problem, not the quality of the approximate zero.nEXERCISESUnless otherwise indicated, use ABSERR = 10-8 andRELERR = 10-6 for the computations with Zero.4.10 Use the code Zero with an initial bracket of [0, l] tofind the roots of the equation F(x) = 0, where F(x) isgiven by each of the following.(a) cos 2x(b) (x - 0.3)(x - 0.6)(c) (x - 0.3)(x - 0.6)(x - 0.9)(d) (x + 1)(x - 0.8)7(e) (x + l)[x7 - 7(0.8)x6 + 21(O.8)2x 5- 35(0.8)3x4 + 35(0.8)4x3 - 21(0.8)5x 2+ 7(0.8)6x - (0.8)7](f) l / cos 2x(g)(h) (x - 3)(x + 1)Print out the approximate roots, FLAG, the number offunction evaluations required, and the residual. Discuss the results. [Sketches of F(x) will be helpful.]4.11 A wire weighs 0.518 lb/ft and is suspended betweentwo towers of equal height (at the same level) and500 ft apart. If the sag in the wire is 50 ft, findthe maximum tension T in the wire.
The appropriateequations to be solved are156CHAPTER 4ROOTS OF NONLINEAR EQUATIONS4.12 For turbulent flow of fluid in a smooth pipe, the equa- 4.17 In trying to solve the equations of radiative transfer intionsemi-infinite atmospheres, one encounters the nonlinear equationmodels the relationship between the friction factor cfand the Reynold’s number Re. Compute cf for Re =104, 105, 106. Solve for all values of the Reynold’snumber in the same run. Do this by communicatingthe parameter Re to the function subprogram using labeled COMMON in FORTRAN or a globally definedvariable in C or C++.where the number ω0, 0 < ω0 < 1, is called an albedo.Show that for fixed ω0, if k is a root, so is -k, and thatthere is a unique value of k with 0 < k < 1 satisfyingthe equation.
For ω0 = 0.25, 0.50, and 0.75, find thecorresponding positive k values. Make some sketchesto help you locate the roots.4.13 In [12] the study of neutron transport in a rod leads toa transcendental equation that has roots related to the 4.18 Exercise 5.27 concerns a temperature distributionproblem where it is necessary to find positive rootscritical lengths. For a rod of length the equation isofMake a rough sketch of the two functions cot(&) and(x2 - l)/(2x) to get an idea of where they intersect toyield roots. For= 1, determine the three smallestpositive roots.4.14 An equation determining the critical load for columnswith batten plates is derived in [9, p.
151]. Suitablevalues of the physical parameters for experiments performed by Timoshenko lead to the problemand the smallest positive root is desired. Make a roughsketch of the function to get an idea of where the rootis. Scale to avoid difficulties with poles and the apparent singularity at 0, and then compute the root.4.15 An equation for the temperature T at which otoluidine has a vapor pressure of 500 mm Hg is foundin [8, p. 424]. In degrees absolute T satisfiesIt is not obvious where the roots are, but a little analysis will help you locate them.
Using brackets fromyour analysis, compute all the roots.4.16 The geometrical concentration factor C in a certain solar energy collection model [10, p. 33] satisfieswhere J0(x) and J1(x) are zeroth and first order Besselfunctions of the first kind. Compute the three smallestpositive roots.4.19 An analysis of the Schrodinger equation for a particle of mass m in a rectangular potential well leads todiscrete sets of values of the total energy E that aresolutions of a pair of transcendental equations. One ofthese equations iswhereis Planck’s constant. Find the value of E that satisfiesthis equation. Use the following data, which correspond to a simplified model of the hydrogen atom:On some machines it may be necessary to scale someof the variables to avoid underflow.
Also, be carefulwith your choice of ABSERR if you want an accurateanswer.4.20 The following problem concerns the cooling of asphere. Suppose the sphere is of radius a and is initially at a temperature V. It cools by Newton’s lawRescale the problem to avoid poles. Find the smallof cooling with thermal conductivity k, thalpance ε,est positive root A if h = 300, C = 1200, f = 0.8, andand diffusivity h2 after being suddenly placed in airD = 14.4.4 CONDITION, LIMITING PRECISION, AND MULTIPLE ROOTSat 0°C.
It can be shown that the temperature θ(r,t) attime t > 0 and radius r isHere the γn are the (positive) roots ofandFor a steel sphere cooling in air at 0°C, suppose theinitial temperature is V = 100°C and the radius isa = 0.30 meters. Appropriate physical constants areh2 = 1.73 × 10-5, ε = 20, and k = 60. Find the threesmallest values of γna and use them to compute A1,A2, and A3. Approximate the temperature at r = 0.25for t = 10 k seconds, k = 2, 3, 4, 5.4.21 When solving f(x) = 0, the subroutine Zero requiresyou to input B and C such that f(B)f(C) < 0. Oftenit is not obvious what to use for B and C, so manyroutines that are similar to Zero begin with a search4.4157for B and C that provide a change of sign.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
