726 Chapter 16.Integration of Ordinary Differential Equations not at all infinitesimal-distance.That single step is a grand leap consisting of many (e.g.,dozens to hundreds)substeps of modified midpoint method,which are then extrapolated to zero stepsize. The sequence of separate attempts to cross the interval H is made with increasing values of n,the number of substeps.Bulirsch and Stoer originally proposed the sequence n=2,4,6,8,12,16,24,32,48,64,96,,[nj=2mj-2, (16.4.1) More recent work by Deufthard [2.3]suggests that the sequence 81 n=2,4,6,8,10,12,14,,[nj=2,. (16.4.2) is usually more efficient.For each step,we do not know in advance how far up this sequence we will go.After each successive n is tried,a polynomial extrapolation ICAL is attempted.That extrapolation gives both extrapolated values and error estimates. If the errors are not satisfactory,we go higher in n.If they are satisfactory,we go on to the next step and begin anew with n =2. RECIPES Of course there must be some upper limit,beyond which we conclude that there is some obstacle in our path in the interval H,so that we must reduce H rather than 9 just subdivide it more finely.In the implementations below,the maximum number of n's to be tried is called KMAXX.For reasons described below we usually take this equal to 8;the 8th value of the sequence (16.4.2)is 16,so this is the maximum 9 number of subdivisions of H that we allow. We enforce error control,as in the Runge-Kutta method,by monitoring internal s爱色日 consistency,and adapting stepsize to match a prescribed bound on the local truncation error.Each new result from the sequence of modified midpoint integrations allows a tableau like that in 83.1 to be extended by one additional set of diagonals.The size of 61 the new correction added at each stage is taken as the (conservative)error estimate How should we use this error estimate to adjust the stepsize?The best strategy now known is due to Deufthard [2,31.For completeness we describe it here: Suppose the absolute value of the error estimate returned from the kth column (and hence the k+1st row)of the extrapolation tableau is e+1.Error control is enforced by requiring Numerica 10621 Ek+1,k <E (16.4.3) 43106 as the criterion for accepting the current step,where e is the required tolerance.For the even sequence (16.4.2)the order of the method is 2+1: k+1,kH2k+1 (16.4.4) Thus a simple estimate ofa new stepsize Hk to obtain convergence in a fixed column k would be 1/(2k+1) Hk=H (16.4.5) Ek+1.k 邑 Which column k should we aim to achieve convergence in?Let's compare the work required for different k.Suppose Ak is the work to obtain row k of the extrapolation tableau, so Ak+1 is the work to obtain column k.We will assume the work is dominated by the cost of evaluating the functions defining the right-hand sides of the differential equations.For nk subdivisions in H,the number of function evaluations can be found from the recurrence A1=n1+1 (16.4.6) Ak+1=Ak+nk+1726 Chapter 16. Integration of Ordinary Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). — not at all infinitesimal — distance. That single step is a grand leap consisting of many (e.g., dozens to hundreds) substeps of modified midpoint method, which are then extrapolated to zero stepsize. The sequence of separate attempts to cross the interval H is made with increasing values of n, the number of substeps. Bulirsch and Stoer originally proposed the sequence n = 2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96,..., [nj = 2nj−2],... (16.4.1) More recent work by Deuflhard [2,3] suggests that the sequence n = 2, 4, 6, 8, 10, 12, 14,..., [nj = 2j],... (16.4.2) is usually more efficient. For each step, we do not know in advance how far up this sequence we will go. After each successive n is tried, a polynomial extrapolation is attempted. That extrapolation gives both extrapolated values and error estimates. If the errors are not satisfactory, we go higher in n. If they are satisfactory, we go on to the next step and begin anew with n = 2. Of course there must be some upper limit, beyond which we conclude that there is some obstacle in our path in the interval H, so that we must reduce H rather than just subdivide it more finely. In the implementations below, the maximum number of n’s to be tried is called KMAXX. For reasons described below we usually take this equal to 8; the 8th value of the sequence (16.4.2) is 16, so this is the maximum number of subdivisions of H that we allow. We enforce error control, as in the Runge-Kutta method, by monitoring internal consistency, and adapting stepsize to match a prescribed bound on the local truncation error. Each new result from the sequence of modified midpoint integrations allows a tableau like that in §3.1 to be extended by one additional set of diagonals. The size of the new correction added at each stage is taken as the (conservative) error estimate. How should we use this error estimate to adjust the stepsize? The best strategy now known is due to Deuflhard [2,3]. For completeness we describe it here: Suppose the absolute value of the error estimate returned from the kth column (and hence the k + 1st row) of the extrapolation tableau is k+1,k. Error control is enforced by requiring k+1,k < (16.4.3) as the criterion for accepting the current step, where is the required tolerance. For the even sequence (16.4.2) the order of the method is 2k + 1: k+1,k ∼ H2k+1 (16.4.4) Thus a simple estimate of a new stepsize Hk to obtain convergence in a fixed column k would be Hk = H
k+1,k 1/(2k+1) (16.4.5) Which column k should we aim to achieve convergence in? Let’s compare the work required for different k. Suppose Ak is the work to obtain row k of the extrapolation tableau, so Ak+1 is the work to obtain column k. We will assume the work is dominated by the cost of evaluating the functions defining the right-hand sides of the differential equations. For nk subdivisions in H, the number of function evaluations can be found from the recurrence A1 = n1 + 1 Ak+1 = Ak + nk+1 (16.4.6)