Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 42
Текст из файла (страница 42)
To illustrate the approach, we refine the Stirling approximation for thecomputation of a large factorial, a simple example that allows the easy computation ofa reference value.REFERENCES207As pointed out in Chapter 1, if we are to evaluate the gamma function for a largeargument, we must scale it to avoid overflow. The Stirling approximation derivedabove suggests that we approximatewhich is both simple and well scaled.
If, for example, we take s = 201, we can obtaina reference value fromThe Stirling approximation iswhich is about 0.17680393 and would beadequate for many purposes. An attempt to approximate the integral simply by truncating the semi-infinite interval to [0, 100] and calling Adapt with ABSERR = 10-8and RELERR = 10-8 results in a value of 0 at a cost of seven function evaluations.If you should find that Adapt requires only seven function evaluations, the smallestnumber it makes, you should give some thought to your problem. Either the integralis very easy or the scale of the problem has been missed entirely. For this integrandx 0 = 1 andThe derivation of the Laplace formula suggeststhat the bulk of the integral is accounted for by the interval [ x0 - 3σ, x0 + 3σ], so wecompute it first and then add on approximations to the integrals over subintervals oflength σ on each side.
The table shows the approximations computed in this way forthe interval [1 - kσ, 1 + kσ] along with the cumulative cost in function evaluations.The last interval here was about [0.5,1.5], making the point that only a small portionof the interval of integration is important. Notice that the code is doing the minimalnumber of function evaluations (14 for the two pieces) for the larger k. The quadratures are very cheap then because the integrand is small and there is an absolute errortolerance.REFERENCES1.
D. Amos, “Evaluation of some cumulative distribution functions by numerical quadrature,”SIAM Review, 20 (1978), pp. 778-800.2. H. Bateman, Partial Differential Equations of Mathematical Physics, Cambridge UniversityPress, London, 1964.3. B. Camahan, H. Luther, and J. Wilkes, Applied Numerical Methods, Wiley, New York, 1969.208CHAPTER 5NUMERICAL INTEGRATION4. P. Davis and P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, NewYork, 1984.5. I. Percival and D. Richards, Introduction to Dynamics, Cambridge University Press, Cambridge,Mass., 1982.6.
F. Fritsch, D. Kahaner, and J. Lyness, “Double integration using one-dimensional adaptivequadrature routines: A software interface problem,” ACM Trans. on Math. Software, 7 (1981),pp. 46-75.7. Handbook of Mathematical Functions, M. Abramowitz and I. Stegun, eds., Dover, Mineola,N.Y., 1964.8. E. Issacson and H. Keller, Analysis of Numerical Methods, Dover, Mineola, N.Y., 1994.9.
V. Kourganoff, Basic Methods in Transfer Problems, Relative Equilibrium and Newton Dimsion, Dover, New York, 1963.10. A. Kronrod, Nodes and Weights of Quadrature Formulas, (Trans.) Consultants Bureau, NewYork, 1965.11. J. Lyness, “When not to use an automatic quadrature routine,” SIAM Review, 25 (1983), pp. 6387.12. R. Piessens, E. de Doncker-Kapenga, C. Überhuber, and D. Kahaner, QUADPACK: SubroutinePackage for Automatic Integration, Springer-Verlag, New York, 1983.13.
A. Stroud and D. Secrest, Gaussian Quadrature Formulas, Prentice Hall, Englewood Cliffs,N.J., 1966.EXERCISESUnless otherwise indicated, use ABSERR = 10-8 andRELERR = 10-6 for the computing exercises.5.25 The potential inside the unit circle due to a specifiedpotential f(θ) on the boundary is given by Poisson’sintegral:There is difficulty evaluating this integral assince for θ´ = θ,and argue that it should have somewhat better numerical properties. Explore this by evaluating φ(r,θ) for rapproaching 1 with f(θ) = sinθ. The analytical solution is then just ϕ(r, θ) = srinθ.5.26 The potential in a conducting strip of width b withpotential zero on the bottom edge and a specified potential F(x) on the upper edge isThis is not too severe because the term is large onlyif r is very close to 1, but in principle there should beno difficulty at all since asBateman [2, pp.
239-241]). Realizing thatwhere β = πy/b. Suppose that an experimenter applies the potential F(x) = 1 for |x| < 0.99 and F(x) =exp[- 100( |x| - 0.99)] for |x| > 0.99. When b = π,compute and plot the potential along the middle of thestrip, ψ(x, π/2).Realizing thatderive the formbound the effect on ϕ(x,y) for yinfinite interval by a finite one:0 of replacing theMiscellaneous Exercises 209For a suitable choice of z, use this instead of the infi- 5.28 If 6what is C? To answer this,nite interval. Show for the given F(x) that ϕ(x,y) =apply Zero with modest tolerances to compute the rootϕ(-x,y), so only a plot for x > 0 is necessary.C of f(C) = 6Evaluate the integral appearing in the function with Adapt and more5.27 This exercise is representative of a great many comstringent tolerances. A little analysis will provide anputations arising in the use of the classical separaappropriate bracket to use for the root solver.
Showtion of variables technique for solving field problems.analytically that the root is simple. This is an examTypically, one must compute many roots of nonlinearple of computing the roots of a function that is relaequations and integrals. The temperature distributiontively expensive to evaluate. How many evaluationsin a cylinder of radius a and height b with the bottomdid Zero need?held at a temperature zero, the top at a temperaturef(r), and the side dissipating heat according to New- 5.29 Reference [3] contains the following problem.
Theton’s law of cooling, can be represented by a series.length L of a carbon dioxide heat exchanger with inIf the thermal conductivity of the cylinder is k and theput temperature T1 and output T0 satisfiesthalpance is ε, then the temperature ϕ(r,z) isThe numbers qn are the positive roots of the equationwhere m = 22.5, D = 0.495, Ts = 550 (temperature ofthe CO2),where the function J0(x) and J1(x) are Bessel functions of the first kind of orders zero and one, respectively. The coefficients An are given byFor T1 = 60, and L = 10, use Adapt and Zero to comThe roots qn depend only on the geometry and the mapute T0.terial. Once they have been computed, one can consider virtually any temperature distribution f(r) by 5.30 The code Adapt uses a queue to hold the subintervals whose quadrature errors are deemed too large tocomputing the quantities An. For k/ε a = 2, we givepass the tolerance test. Do you think there would bethe problem of solving for qna for n = 1, 2, 3 in Exany difference in performance if a stack were used inercise 4.18.
If you have not done that problem, thestead? Explain. Modify Adapt so that it uses a stackroots are 0.940770563949734, 3.95937118501253,and test it on some difficult problems to see what acand 7.08638084796766. Then for a = 1, computetually happens.A 1, A2, A3 for f(r) = exp(-r) - exp(-1).Previous HomeCHAPTER 6ORDINARY DIFFERENTIALEQUATIONSHistorically, ordinary differential equations have originated in chemistry, physics, andengineering.
They still do, but nowadays they also originate in medicine, biology,anthropology, and the like. Because differential equations are more complex theoretically than the other topics taken up in this book and because their numerical solutionis less well understood, the art of numerical computing plays an important role in thischapter. There is a subdiscipline of numerical analysis called mathematical softwarethat colors all the work of this book. It is particularly evident in this chapter with itsmore complex task.6.1SOME ELEMENTS OF THE THEORYThe simplest kind of ordinary differential equation problem is to find a function y(x)with a derivative continuous on [a,b] such thaty´(x) = f(x),a < x < b,for a given continuous function f(x).
From elementary calculus it is known that such afunction y(x) exists—the indefinite integral of f(x). However, if y(x) satisfies the differential equation, so does y(x) + c for any constant c. To specify a particular solution,some more information about y(x) is required. The most common kind of additionalinformation supplied is the value A of y(x) at the initial point a. Thenis the unique solution to the initial value problem consisting ofsatisfied by y(x) and the initial value of y(x).The general first order ordinary differential equation has fas x.
It is assumed that f(x,y) is continuous for a < x < b anda function of x with a continuous derivative for a < x < b thaty´(x) = f(x,y(x))210the differential equationdepending on y as wellall y. A solution y(x) issatisfies the equation(6.1)Next6.1 SOME ELEMENTS OF THE THEORY211for each x in the interval [a,b]. Typically the solution desired is specified by its valueat the initial point of the interval:y(a) = A.(6.2)Equation (6.2) is called an initial condition, and the combination of (6.1) and (6.2) iscalled an initial value problem for an ordinary differential equation.In elementary treatments of differential equations, the initial value problem has aunique solution that exists throughout the interval of interest and that can be obtainedby analytical means (more familiarly called a “trick”).
However, for most problemsthat are not contrived, an analytical solution is impossible to obtain or is less satisfactory than a numerical solution. Matters are also complicated by the fact that solutionscan fail to exist over the desired interval of interest. Problems with solutions that“blow up” place a special burden on a numerical procedure, although we might wellexpect a general-purpose code to compute such solutions until overflow occurs. Problems that have more than one solution are especially troublesome. Difficulties withexistence and uniqueness will be excluded at the level of the theory to be developedhere.
A simple condition that guarantees that these difficulties will not occur can beformulated in terms of how f(x,y) depends on y.The function f(x,y) satisfies a Lipschitz condition in y if for all x in the interval[a,b] and for all u, v,with L a constant, hereafter called a Lipschitz constant. The inequality assumes a morefamiliar form if f has a continuous partial derivative in its second variable, for thenfor some w between u and v. Ifis bounded in magnitude for all arguments, thenf satisfies a Lipschitz condition and any constant L such thatfor all x in [a,b] and all w is a Lipschitz constant.
If the partial derivative is notbounded, it is not hard to show that the inequality (6.3) cannot hold for all u, v and allx in [a,b], so f does not satisfy a Lipschitz condition.Example 6.1. The function f(x,y) = x2 cos2 y + y sin2 x defined for |x| < 1 and ally satisfies a Lipschitz condition with constant L = 3. To see this, differentiate withrespect to y to getand so for the range of x allowedn212CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSExample 4.2. The function f(x,y)does not satisfy a Lipschitz conditionbecause it has a continuous partial derivative for y > 0 that is not bounded asAn important special case of (6.1) is that of a linear differential equation, an equation of the form f(x,y) = g(x)y + h(x). The function f(x,y) being continuous in (x,y)is then equivalent to g(x) and h(x) being continuous in x. Becauseand because a continuous function g(x) is bounded in magnitude on any finite interval[a,b], a linear equation is Lipschitzian in nearly all cases of practical interest.Example 6.3.Dawson’s integral is the functionYou should verify that it is a solution of the initial value problem for the linear differential equationy´ = l - 2x yy(0) = 0.On the interval [0,b] for any b 0, the function f(x,y) = 1 - 2 xy is continuous andLipschitzian with Lipschitz constant L = 2|b|.nSufficient conditions for existence and uniqueness can now be stated formally.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
