16.4 Richardson Extrapolation and the Bulirsch-Stoer Method 727 The work per unit step to get column k is A/Hk,which we nondimensionalize with a factor of H and write as WE =Ak+LH (16.4.7) Hk =Ak+(+达)e+) (16.4.8) The quantities W can be calculated during the integration.The optimal column index g is then defined by Wg=kmi血，Wa (16.4.9 .....f where ky is the final column,in which the error criterion(16.4.3)was satisfied.The g determined from(16.4.9)defines the stepsize H to be used as the next basic stepsize,so that we can expect to get convergence in the optimal column g. Two important refinements have to be made to the strategy outlined so far: ICAL ·If the current H is“too small,”then kr will be“too small,.”and so g remains "too small."It may be desirable to increase H and aim for convergence in a 3 column g kf. If the current H is"too big,"we may not converge at all on the current step and we RECIPES will have to decrease H.We would like to detect this by monitoring the quantities ek+1.for each k so we can stop the current step as soon as possible. Deufhard's prescription for dealing with these two problems uses ideas from communi- cation theory to determine the "average expected convergence behavior"of the extrapolation. 传 令 Press. His model produces certain correction factors a(k,g)by which Hk is to be multiplied to try to get convergence in column q.The factors a(k,q)depend only on e and the sequence {n and so can be computed once during initialization: ak,）=e可 for k<g (16.4.10) with a(g,q)=1. 6 Now to handle the first problem,suppose convergence occurs in column g=ky.Then g rather than taking H for the next step,we might aim to increase the stepsize to get convergence in column q+1.Since we don't have H+available from the computation,we estimate it as Hg+1=Hgo(g,q+1) (16.4.11) By equation (16.4.7)this replacement is efficient,i.e.,reduces the work per unit step,if ◇ Ag+1> Aq+2 Numerical Ha (16.4.12) -43126 or Ag+ia(q;q+1)>Ag+2 (16.4.13) During initialization,this inequality can be checked for g=1,2,...to determine kmax,the largest allowed column.Then when (16.4.12)is satisfied it will always be efficient to use H+1.(In practice we limit max to 8 even when e is very small as there is very little further gain in efficiency whereas roundoff can become a problem.) The problem of stepsize reduction is handled by computing stepsize estimates ik≡Hka(k,q), k=1,..,q-1 (16.4.14) during the current step.The 's are estimates of the stepsize to get convergence in the optimal column g.If any H is "too small,"we abandon the current step and restart using H.The criterion of being“too small'”is taken to be Hka(k,q+1<H (16.4.15) The a's satisfy a(k,q+1)>a(k,q).16.4 Richardson Extrapolation and the Bulirsch-Stoer Method 727 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). The work per unit step to get column k is Ak+1/Hk, which we nondimensionalize with a factor of H and write as Wk = Ak+1 Hk H (16.4.7) = Ak+1 k+1,k 1/(2k+1) (16.4.8) The quantities Wk can be calculated during the integration. The optimal column index q is then defined by Wq = min k=1,...,kf Wk (16.4.9) where kf is the final column, in which the error criterion (16.4.3) was satisfied. The q determined from (16.4.9) defines the stepsize Hq to be used as the next basic stepsize, so that we can expect to get convergence in the optimal column q. Two important refinements have to be made to the strategy outlined so far: • If the current H is “too small,” then kf will be “too small,” and so q remains “too small.” It may be desirable to increase H and aim for convergence in a column q>kf . • If the current H is “too big,” we may not converge at all on the current step and we will have to decrease H. We would like to detect this by monitoring the quantities k+1,k for each k so we can stop the current step as soon as possible. Deuflhard’s prescription for dealing with these two problems uses ideas from communi￾cation theory to determine the “average expected convergence behavior” of the extrapolation. His model produces certain correction factors α(k, q) by which Hk is to be multiplied to try to get convergence in column q. The factors α(k, q) depend only on and the sequence {ni} and so can be computed once during initialization: α(k, q) = Ak+1−Aq+1 (2k+1)(Aq+1−A1+1) for k<q (16.4.10) with α(q, q)=1. Now to handle the first problem, suppose convergence occurs in column q = kf . Then rather than taking Hq for the next step, we might aim to increase the stepsize to get convergence in column q + 1. Since we don’t have Hq+1 available from the computation, we estimate it as Hq+1 = Hqα(q, q + 1) (16.4.11) By equation (16.4.7) this replacement is efficient, i.e., reduces the work per unit step, if Aq+1 Hq > Aq+2 Hq+1 (16.4.12) or Aq+1α(q, q + 1) > Aq+2 (16.4.13) During initialization, this inequality can be checked for q = 1, 2,... to determine kmax, the largest allowed column. Then when (16.4.12) is satisfied it will always be efficient to use Hq+1. (In practice we limit kmax to 8 even when is very small as there is very little further gain in efficiency whereas roundoff can become a problem.) The problem of stepsize reduction is handled by computing stepsize estimates H¯k ≡ Hkα(k, q), k = 1,...,q − 1 (16.4.14) during the current step. The H¯ ’s are estimates of the stepsize to get convergence in the optimal column q. If any H¯k is “too small,” we abandon the current step and restart using H¯k. The criterion of being “too small” is taken to be Hkα(k, q + 1) < H (16.4.15) The α’s satisfy α(k, q + 1) > α(k, q).