710 Chapter 16.Integration of Ordinary Differential Equations CITED REFERENCES AND FURTHER READING: Gear,C.W.1971,Numerical Initial Value Problems in Ordinary Differential Equations(Englewood Cliffs,NJ:Prentice-Hall). Acton,F.S.1970.Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 5. Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Chapter 7. Lambert,J.1973,Computational Methods in Ordinary Differential Equations(New York:Wiley). Lapidus,L.,and Seinfeld,J.1971,Numerical Solution of Ordinary Differential Eguations (New York:Academic Press). 16.1 Runge-Kutta Method The formula for the Euler method is ⊙ Un+1 =Un hf(In;yn) (16.1.1) which advances a solution fromn to=n+h.The formula is unsymmetrical: % Press. It advances the solution through an interval h,but uses derivative information only at the beginning of that interval(see Figure 16.1.1).That means(and you can verify by expansion in power series)that the step's error is only one power of h smaller ◆分 than the correction,i.e O(h2)added to (16.1.1). There are several reasons that Euler's method is not recommended for practical use,among them,(i)the method is not very accurate when compared to other, 61 fancier,methods run at the equivalent stepsize,and(ii)neither is it very stable (see $16.6 below). Consider,however,the use of a step like (16.1.1)to take a"trial"step to the midpoint of the interval.Then use the value of both x and y at that midpoint to compute the "real"step across the whole interval.Figure 16.1.2 illustrates the idea.In equations, k1=hf(In,Un) 、6心 9 Numerica 10621 43106 k2 hf (In+5h,yn +k1) (16.1.2) (outside n+1=n+k2+O(h3) North As indicated in the error term,this symmetrization cancels out the first-order error term,making the method second order.[A method is conventionally called nth order if its error term is O(hn+1).]In fact,(16.1.2)is called the second-order Runge-Kutta or midpoint method. We needn't stop there.There are many ways to evaluate the right-hand side f(,y)that all agree to first order,but that have different coefficients of higher-order error terms.Adding up the right combination of these,we can eliminate the error terms order by order.That is the basic idea of the Runge-Kutta method.Abramowitz and Stegun [1],and Gear [2],give various specific formulas that derive from this basic
710 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 machinereadable 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). CITED REFERENCES AND FURTHER READING: Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (Englewood Cliffs, NJ: Prentice-Hall). Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 5. Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 7. Lambert, J. 1973, Computational Methods in Ordinary Differential Equations (New York: Wiley). Lapidus, L., and Seinfeld, J. 1971, Numerical Solution of Ordinary Differential Equations (New York: Academic Press). 16.1 Runge-Kutta Method The formula for the Euler method is yn+1 = yn + hf(xn, yn) (16.1.1) which advances a solution from xn to xn+1 ≡ xn+h. The formula is unsymmetrical: It advances the solution through an interval h, but uses derivative information only at the beginning of that interval (see Figure 16.1.1). That means (and you can verify by expansion in power series) that the step’s error is only one power of h smaller than the correction, i.e O(h2) added to (16.1.1). There are several reasons that Euler’s method is not recommended for practical use, among them, (i) the method is not very accurate when compared to other, fancier, methods run at the equivalent stepsize, and (ii) neither is it very stable (see §16.6 below). Consider, however, the use of a step like (16.1.1) to take a “trial” step to the midpoint of the interval. Then use the value of both x and y at that midpoint to compute the “real” step across the whole interval. Figure 16.1.2 illustrates the idea. In equations, k1 = hf(xn, yn) k2 = hf xn + 1 2h, yn + 1 2 k1 yn+1 = yn + k2 + O(h3) (16.1.2) As indicated in the error term, this symmetrization cancels out the first-order error term, making the method second order. [A method is conventionally called nth order if its error term is O(hn+1).] In fact, (16.1.2) is called the second-order Runge-Kutta or midpoint method. We needn’t stop there. There are many ways to evaluate the right-hand side f(x, y) that all agree to first order, but that have different coefficients of higher-order error terms. Adding up the right combination of these, we can eliminate the error terms order by order. That is the basic idea of the Runge-Kutta method. Abramowitz and Stegun [1], and Gear [2], give various specific formulas that derive from this basic
16.1 Runge-Kutta Method 711 Mx) ② ① X2 x 83g Figure 16.1.1.Euler's method.In this simplest (and least accurate)method for integrating an ODE, granted for 19881992 the derivative at the starting point of each interval is extrapolated to find the next function value.The method has first-order accuracy. 1-800 y(x) to any users to Cambridge from NUMERICAL RECIPES I 4 ① t-- University Press. THE 3 5 (North America 州bMe se make one paper ART ictly proh Programs X2 X3 Figure 16.1.2.Midpoint method.Second-order accuracy is obtained by using the initial derivative at each step to find a point halfway across the interval,then using the midpoint derivative across the full width of the interval.In the figure,filled dots represent final function values,while open dots represent to dir function values that are discarded once their derivatives have been calculated and used. idea.By far the most often used is the classical fourth-order Runge-Kutta formula, 1788-1982 OF SCIENTIFIC COMPUTING(ISBN which has a certain sleekness of organization about it: k1 hf(In:Un) Numerical Recipes 10-521 h k2 hf(n+2:Un+2 43108 h k2 ks hf(n+2:Un2 (outside 膜 ka hf(n +h,yn +k3) Software. 2大3 KA Un+1 yn+ 6 +3 +O(h5) (16.1.3) ying of The fourth-order Runge-Kutta method requires four evaluations of the right- hand side per step h(see Figure 16.1.3).This will be superior to the midpoint method(16.1.2)ifat least twice as large a step is possible with(16.1.3)for the same accuracy.Is that so?The answer is:often,perhaps even usually,but surely not always!This takes us back to a central theme,namely that high order does not always mean high accuracy.The statement"fourth-order Runge-Kutta is generally superior to second-order"is a true one,but you should recognize it as a statement about the
16.1 Runge-Kutta Method 711 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 machinereadable 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). y(x) 1 2 x1 x2 x3 x Figure 16.1.1. Euler’s method. In this simplest (and least accurate) method for integrating an ODE, the derivative at the starting point of each interval is extrapolated to find the next function value. The method has first-order accuracy. y(x) 1 2 x1 x2 x3 x 3 4 5 Figure 16.1.2. Midpoint method. Second-order accuracy is obtained by using the initial derivative at each step to find a point halfway across the interval, then using the midpoint derivative across the full width of the interval. In the figure, filled dots represent final function values, while open dots represent function values that are discarded once their derivatives have been calculated and used. idea. By far the most often used is the classical fourth-order Runge-Kutta formula, which has a certain sleekness of organization about it: k1 = hf(xn, yn) k2 = hf(xn + h 2 , yn + k1 2 ) k3 = hf(xn + h 2 , yn + k2 2 ) k4 = hf(xn + h, yn + k3) yn+1 = yn + k1 6 + k2 3 + k3 3 + k4 6 + O(h5) (16.1.3) The fourth-order Runge-Kutta method requires four evaluations of the righthand side per step h (see Figure 16.1.3). This will be superior to the midpoint method (16.1.2) if at least twice as large a step is possible with (16.1.3) for the same accuracy. Is that so? The answer is: often, perhaps even usually, but surely not always! This takes us back to a central theme, namely that high order does not always mean high accuracy. The statement “fourth-order Runge-Kutta is generally superior to second-order” is a true one, but you should recognize it as a statement about the
712 Chapter 16.Integration of Ordinary Differential Equations http://www.nr. Figure 16.1.3.Fourth-order Runge-Kutta method.In each step the derivative is evaluated four times: once at the initial point,twice at trial midpoints,and once at a trial endpoint.From these derivatives the 8 final function value (shown as a filled dot)is calculated.(See text for details.) 鱼 18881892 contemporary practice of science rather than as a statement about strict mathematics. 100 That is,it reflects the nature of the problems that contemporary scientists like to solve. For many scientific users,fourth-order Runge-Kutta is not just the first word on ODE integrators,but the last word as well.In fact,you can get pretty far on this old workhorse,especially if you combine it with an adaptive stepsize algorithm.Keep 豆是 令 in mind,however,that the old workhorse's last trip may well be to take you to the poorhouse:Bulirsch-Stoer or predictor-corrector methods can be very much more efficient for problems where very high accuracy is a requirement.Those methods are the high-strung racehorses.Runge-Kutta is for ploughing the fields.However. even the old workhorse is more nimble with new horseshoes.In $16.2 we will give 9 a modern implementation of a Runge-Kutta method that is quite competitive as long as very high accuracy is not required.An excellent discussion of the pitfalls in constructing a good Runge-Kutta code is given in [3] Here is the routine for carrying out one classical Runge-Kutta step on a set of n differential equations.You input the values of the independent variables,and you get out new values which are stepped by a stepsize h(which can be positive or negative).You will notice that the routine requires you to supply not only function 是家 derivs for calculating the right-hand side,but also values of the derivatives at the starting point.Why not let the routine call derivs for this first value?The answer will become clear only in the next section.but in brief is this:This call may not be 10-621 your only one with these starting conditions.You may have taken a previous step with too large a stepsize,and this is your replacement.In that case,you do not 43106 want to call derivs unnecessarily at the start.Note that the routine that follows has,therefore,only three calls to derivs. (outside Software. #include "nrutil.h" North 6 void rk4(float y],float dydx[],int n,float x,float h,float yout[], void (*derivs)(float,float []float []) Given values for the variables y[1..n]and their derivatives dydx [1..n]known at x,use the B fourth-order Runge-Kutta method to advance the solution over an interval h and return the incremented variables as yout[1..n],which need not be a distinct array from y.The user supplies the routine derivs(x,y,dydx),which returns derivatives dydx at x. int i; float xh,hh,h6,*dym,*dyt,*yt; dym=vector(1,n); dyt=vector(1,n);
712 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 machinereadable 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). 1 2 3 4 yn + 1 yn Figure 16.1.3. Fourth-order Runge-Kutta method. In each step the derivative is evaluated four times: once at the initial point, twice at trial midpoints, and once at a trial endpoint. From these derivatives the final function value (shown as a filled dot) is calculated. (See text for details.) contemporary practice of science rather than as a statement about strict mathematics. That is, it reflects the nature of the problems that contemporary scientists like to solve. For many scientific users, fourth-order Runge-Kutta is not just the first word on ODE integrators, but the last word as well. In fact, you can get pretty far on this old workhorse, especially if you combine it with an adaptive stepsize algorithm. Keep in mind, however, that the old workhorse’s last trip may well be to take you to the poorhouse: Bulirsch-Stoer or predictor-corrector methods can be very much more efficient for problems where very high accuracy is a requirement. Those methods are the high-strung racehorses. Runge-Kutta is for ploughing the fields. However, even the old workhorse is more nimble with new horseshoes. In §16.2 we will give a modern implementation of a Runge-Kutta method that is quite competitive as long as very high accuracy is not required. An excellent discussion of the pitfalls in constructing a good Runge-Kutta code is given in [3]. Here is the routine for carrying out one classical Runge-Kutta step on a set of n differential equations. You input the values of the independent variables, and you get out new values which are stepped by a stepsize h (which can be positive or negative). You will notice that the routine requires you to supply not only function derivs for calculating the right-hand side, but also values of the derivatives at the starting point. Why not let the routine call derivs for this first value? The answer will become clear only in the next section, but in brief is this: This call may not be your only one with these starting conditions. You may have taken a previous step with too large a stepsize, and this is your replacement. In that case, you do not want to call derivs unnecessarily at the start. Note that the routine that follows has, therefore, only three calls to derivs. #include "nrutil.h" void rk4(float y[], float dydx[], int n, float x, float h, float yout[], void (*derivs)(float, float [], float [])) Given values for the variables y[1..n] and their derivatives dydx[1..n] known at x, use the fourth-order Runge-Kutta method to advance the solution over an interval h and return the incremented variables as yout[1..n], which need not be a distinct array from y. The user supplies the routine derivs(x,y,dydx), which returns derivatives dydx at x. { int i; float xh,hh,h6,*dym,*dyt,*yt; dym=vector(1,n); dyt=vector(1,n);
16.1 Runge-Kutta Method 713 yt=vector(1,n); hh=h*0.5: h6=h/6.0; xh=x+hh; for (i=1;i<=n;i++)yt [i]=y[i]+hh*dydx[i]; First step. (*derivs)(xh,yt,dyt); Second step. for (i=1;i<=n;i++)yt[i]-y[i]+hh*dyt [i]; (*derivs)(xh,yt,dym); Third step for(i=1:1<=n;i++)[ yt[i]=y[i]+h*dym[i]; dym[i]+dyt[i]; http://www. (*derivs)(x+h,yt,dyt); Fourth step. for(i=1;1<=n;i+) Accumulate increments with proper yout [i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]); weights. free_vector(yt,1,n); 鱼 18881892 free_vector(dyt,1,n); free_vector(dym,1,n); 1800 from NUMERICAL RECIPESI The Runge-Kutta method treats every step in a sequence of steps in identical manner.Prior behavior of a solution is not used in its propagation.This is mathematically proper,since any point along the trajectory of an ordinary differential equation can serve as an initial point.The fact that all steps are treated identically also makes it easy to incorporate Runge-Kutta into relatively simple"driver"schemes. Americ computer Press. We consider adaptive stepsize control,discussed in the next section,an essential ART for serious computing.Occasionally,however,you just want to tabulate a function at 9 Progra equally spaced intervals,and without particularly high accuracy.In the most common case,you want to produce a graph of the function.Then all you need may be a simple driver program that goes from an initial to a final f in a specified number of steps.To check accuracy,double the number of steps,repeat the integration,and compare results.This approach surely does not minimize computer time,and it can to dir fail for problems whose nature requires a variable stepsize,but it may well minimize user effort.On small problems,this may be the paramount consideration. Here is such a driver,self-explanatory,which tabulates the integrated functions in the global arrays *x and **y;be sure to allocate memory for them with the routines vector()and matrix(),respectively. #include "nrutil.h" OF SCIENTIFIC COMPUTING(ISBN 19841820 Numerical Recipes 10-621 -43108 float **y,*xx; For communication back to main. (outside void rkdumb(float vstart[],int nvar,float x1,float x2,int nstep, Software. void (*derivs)(float,float []float []) Starting from initial values vstart [1..nvar]known at x1 use fourth-order Runge-Kutta to advance nstep equal increments to x2.The user-supplied routine derivs(x,v,dvdx) evaluates derivatives.Results are stored in the global variables y[1..nvar][1..nstep+1] and xx[1..nstep+1]. void rk4(float y[],float dydx[],int n,float x,float h,float yout[], void(*derivs)(f1oat,f1oat[□，f1oat[□))； int i,k; float x,hi float *v,*vout,*dv; v=vector(1,nvar); vout-vector(1,nvar);
16.1 Runge-Kutta Method 713 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 machinereadable 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). yt=vector(1,n); hh=h*0.5; h6=h/6.0; xh=x+hh; for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i]; First step. (*derivs)(xh,yt,dyt); Second step. for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i]; (*derivs)(xh,yt,dym); Third step. for (i=1;i<=n;i++) { yt[i]=y[i]+h*dym[i]; dym[i] += dyt[i]; } (*derivs)(x+h,yt,dyt); Fourth step. for (i=1;i<=n;i++) Accumulate increments with proper yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]); weights. free_vector(yt,1,n); free_vector(dyt,1,n); free_vector(dym,1,n); } The Runge-Kutta method treats every step in a sequence of steps in identical manner. Prior behavior of a solution is not used in its propagation. This is mathematically proper, since any point along the trajectory of an ordinary differential equation can serve as an initial point. The fact that all steps are treated identically also makes it easy to incorporate Runge-Kutta into relatively simple “driver” schemes. We consider adaptive stepsize control, discussed in the next section, an essential for serious computing. Occasionally, however, you just want to tabulate a function at equally spaced intervals, and without particularly high accuracy. In the most common case, you want to produce a graph of the function. Then all you need may be a simple driver program that goes from an initial xs to a final xf in a specified number of steps. To check accuracy, double the number of steps, repeat the integration, and compare results. This approach surely does not minimize computer time, and it can fail for problems whose nature requires a variable stepsize, but it may well minimize user effort. On small problems, this may be the paramount consideration. Here is such a driver, self-explanatory, which tabulates the integrated functions in the global arrays *x and **y; be sure to allocate memory for them with the routines vector() and matrix(), respectively. #include "nrutil.h" float **y,*xx; For communication back to main. void rkdumb(float vstart[], int nvar, float x1, float x2, int nstep, void (*derivs)(float, float [], float [])) Starting from initial values vstart[1..nvar] known at x1 use fourth-order Runge-Kutta to advance nstep equal increments to x2. The user-supplied routine derivs(x,v,dvdx) evaluates derivatives. Results are stored in the global variables y[1..nvar][1..nstep+1] and xx[1..nstep+1]. { void rk4(float y[], float dydx[], int n, float x, float h, float yout[], void (*derivs)(float, float [], float [])); int i,k; float x,h; float *v,*vout,*dv; v=vector(1,nvar); vout=vector(1,nvar);
714 Chapter 16.Integration of Ordinary Differential Equations dv=vector(1,nvar); for(1=1;1<=nvar;i++){ Load starting values. v[i]-vstart[i]; y[i][1]=v[1]; 2 xx[1]=x1; x=x1; h=(x2-x1)/nstep; for (k=1;k<=nstep;k++){ Take nstep steps. (*derivs)(x,v,dv); rk4(v,dv,nvar,x,h,vout,derivs); if ((float)(x+h)==x)nrerror("Step size too small in routine rkdumb"); X+=h; xx[0k+1]=x; Store intermediate steps. for (i=1;i<=nvar;i++){ v[i]-vout[i]: y[i][k+1]=v[i]; 2 2 令 1-800 NUMERI ICAL free_vector(dv,1,nvar); free_vector(vout,1,nvar); free_vector(v,1,nvar); RECIPES I (North America computer, Press. CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by 9 Programs Dover Publications,New York),825.5.[1] Gear,C.W.1971,Numerical Initial Value Problems in Ordinary Differential Equations(Englewood 、核 OF SCIENTIFIC Cliffs,NJ:Prentice-Hall),Chapter 2.[2] Shampine,L.F.,and Watts,H.A.1977,in Mathematica/Software /l,J.R.Rice,ed.(New York:Aca- 61 demic Press),pp.257-275:1979,Applied Mathematics and Computation,vol.5,pp.93- 121.[3] Rice,J.R.1983,Numerical Methods,Software,and Analysis (New York:McGraw-Hill),89.2. Numerica 10.621 16.2 Adaptive Stepsize Control for Runge-Kutta 431 Recipes A good ODE integrator should exert some adaptive control over its own progress, (outside making frequent changes in its stepsize.Usually the purpose of this adaptive stepsize 首 control is to achieve some predetermined accuracy in the solution with minimum computational effort.Many small steps should tiptoe through treacherous terrain. while a few great strides should speed through smooth uninteresting countryside. The resulting gains in efficiency are not mere tens of percents or factors of two; they can sometimes be factors of ten,a hundred,or more.Sometimes accuracy may be demanded not directly in the solution itself,but in some related conserved quantity that can be monitored. Implementation of adaptive stepsize control requires that the stepping algorithm signal information about its performance,most important,an estimate of its truncation error.In this section we will learn how such information can be obtained.Obviously
714 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 machinereadable 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). dv=vector(1,nvar); for (i=1;i<=nvar;i++) { Load starting values. v[i]=vstart[i]; y[i][1]=v[i]; } xx[1]=x1; x=x1; h=(x2-x1)/nstep; for (k=1;k<=nstep;k++) { Take nstep steps. (*derivs)(x,v,dv); rk4(v,dv,nvar,x,h,vout,derivs); if ((float)(x+h) == x) nrerror("Step size too small in routine rkdumb"); x += h; xx[k+1]=x; Store intermediate steps. for (i=1;i<=nvar;i++) { v[i]=vout[i]; y[i][k+1]=v[i]; } } free_vector(dv,1,nvar); free_vector(vout,1,nvar); free_vector(v,1,nvar); } CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §25.5. [1] Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 2. [2] Shampine, L.F., and Watts, H.A. 1977, in Mathematical Software III, J.R. Rice, ed. (New York: Academic Press), pp. 257–275; 1979, Applied Mathematics and Computation, vol. 5, pp. 93– 121. [3] Rice, J.R. 1983, Numerical Methods, Software, and Analysis (New York: McGraw-Hill), §9.2. 16.2 Adaptive Stepsize Control for Runge-Kutta A good ODE integrator should exert some adaptive control over its own progress, making frequent changes in its stepsize. Usually the purpose of this adaptive stepsize control is to achieve some predetermined accuracy in the solution with minimum computational effort. Many small steps should tiptoe through treacherous terrain, while a few great strides should speed through smooth uninteresting countryside. The resulting gains in efficiency are not mere tens of percents or factors of two; they can sometimes be factors of ten, a hundred, or more. Sometimes accuracy may be demanded not directly in the solution itself, but in some related conserved quantity that can be monitored. Implementation of adaptive stepsize control requires that the stepping algorithm signal information about its performance,most important, an estimate of its truncation error. In this section we will learn how such information can be obtained. Obviously