24. Ordinary Differential Equations¶
PDF pages 1121–1132
24.1.2 O and o¶
Chapter 24 Ordinary Differential Equations Here we’ll cover some of the basic ideas for numerically solving ordinary differential equations, particularly for solving equations of motion for the dynamical evolution of a system. We won’t be aiming for completeness in this vast area, but rather we will just cover some of the important and useful methods. First, however, we need a few mathematical preliminaries that will allow us to quantify how much of an error we’re making when generating a numerical—and hence approximate—solution to an ODE. 24.1 Convergence 24.1.1 Sequences First, let’s start off with the notion of a sequence. A real-valued sequence is a function x : Z+ −\rightarrow R, where Z+ denotes the set of positive integers. The sequence is commonly written using the subscript notation xn. That is, the sequence associates a real number with every positive integer. The generalization of the real-valued sequence to complex values and n-tuples is obvious in the ideas of convergence to follow, so long as an appropriate ‘‘distance,’’ such as the Euclidean metric, is used.
is an N \in Z+ such that if n > N, then |xn −x| < ϵ. In this case, x is the limit of the sequence xn. If the sequence does not converge to a limit, then the sequence is said to diverge, in which case there are still a couple of possibilities. If for any x \in R, there is an N \in Z+ such that if n > N, then xn > x, then the sequence xn is said to diverge to \infty. Similarly, if for the same conditions xn < x, then the sequence xn is said to diverge to −\infty. If the sequence neither converges nor diverges to \pm\infty, it is said to oscillate. As a quick example, let’s prove formally that the sequence xn = 1/n converges to zero. We start by letting ϵ > 0. Choose N to be the smallest integer larger than 1/ϵ so that 1/N < ϵ. If n > N, then
24.1.2 O and o Now supposing that we have a sequence, we can address the question of how the sequence behaves asymp- totically by comparing the two sequences. Thus, let xn and yn be two sequences, not necessarily convergent. Then we write
(24.1) read as ‘‘yn is big Oh of xn’’ or ‘‘yn is of the order of xn,’’ if there is a K > 0 and an N \in Z+, such that if n > N, then |yn/xn| < K. On the other hand, if the sequence (yn/xn) −\rightarrow 0, then yn evidently converges more quickly than xn, and we write
(24.2) This statement is read as ‘‘yn is little oh of xn,’’ or ‘‘yn is dominated by xn.’’ Note that the notation
24.1.3 Convergence and Scaling of Functions¶
Chapter 24. Ordinary Differential Equations 24.1.2.1 Example To illustrate these definitions, consider the sequence xn = 3 + 3n + 5n2.
24.1.3 Convergence and Scaling of Functions These same notions of convergence also apply to limits of functions. In particular, when we say that lim
(24.3) we mean that for every sequence xn −\rightarrow x0, the corresponding sequence f(xn) −\rightarrow y. If it turns out that
than the limit, xn −x = |xn −x|, then the same definition gives the limit from above lim x\rightarrow x+
(24.4) while if we consider only sequences of numbers smaller than the limit, x −xn = |x −xn|, then the same definition gives the limit from below lim x\rightarrow x−
(24.5) assuming these limits exist. The function is then continuous if the limits agree, y+ = y−. The concepts of O and o can similarly be extended to functions of a continuum. In particular, we can write
as x −\rightarrow x0 (24.6)
the condition that each sequence |f(xn)/g(xn)| is eventually bounded, lim n\rightarrow \infty
f(xn) g(xn) < \infty. (24.7) and thus we may equivalently write lim x\rightarrow x0
f(x) g(x) < \infty. (24.8) For the little o, we can similarly write
(24.9)
there is only a sense of ‘‘little oh’’ dominance of functions as the argument increases without bound. 24.1.3.1 Truncation Error The O notation is particularly useful for representing the error in a truncated series expansion, or the truncation error. That is, suppose we write the series expansion of the exponential function near x = 0 as
2! + O(x3). (24.10) Since the remaining terms in the expansion are −x3/3!, x4/4!, and so on, we can compactly represent them by simply noting that all the unwritten terms are O(x3). This gives a rough idea of how the truncation error scales with x. The same idea obviously applies to asymptotic expansions about the point at \infty.
24.2.1 Local and Global Truncation Error¶
24.2 Euler Methods 24.1.3.2 Example¶
24.2 Euler Methods Consider an ordinary differential equation of the form¶
(24.11) that we wish to solve. The Euler method is the simplest way to generate the solution y(t) to the ODE. The idea is to consider the Taylor expansion of the evolved solution
(24.12) Thus, the Euler method consists of making the approximation
(24.13) (Euler method) to obtain the advanced solution y(t + ∆t), since the function f and the present solution y(t) are known. This process is iterated to generate the further advanced solutions y(t + 2∆t), y(t + 3∆t), and so on. Pretty simple, and pretty easy to implement on a computer. This is the simplest example of a finite-difference method, since finite time steps are taken to approximate the continuous solution. In a slightly more compact notation, we may write the Euler method as the recurrence equation
(24.14) (Euler method)
Of course, this same method can be applied to system of coupled ODEs. If we have a set of ODEs of the form
(24.15) then the finite-step update equations are
(24.16) which is a pretty obvious generalization of the scalar case. In this sense, the Euler method may be applied generally to any problem, since any ODE can be decomposed into a system of first-order ODEs (an nth- order ODE decomposes into n coupled, first-order ODEs). However, typically a somewhat fancier method is usually appropriate. 24.2.1 Local and Global Truncation Error Evidently, the error in taking a single step in the Euler method is O(∆t2), as we see from Eq. (24.12). Thus, we say that the local trunction error of the Euler method is O(∆t2). The global truncation error refers to the error in generating the solution over a fixed interval, say from 0 to t in time. This takes N = t/∆t steps, and in the worst case when the local truncation errors add,
24.2.2 Implicit Euler Method and Stiff Equations¶
Chapter 24. Ordinary Differential Equations
the errors are not random, but are predictably related to the form of f(y(t), t) and the solution y(t). In any case, in principle the solution converges to the true one as ∆t −\rightarrow 0. The Euler method is said to be a first-order method, meaning that the local truncation is correct to first order in ∆t, or that the global truncation error is first order in ∆t. 24.2.2 Implicit Euler Method and Stiff Equations An alternative to the expansion in Eq. (24.12) is to expand y(t) about t + ∆t:
(24.17) Thus, an alternative stepping scheme is the implicit Euler method, given by rearranging the expansion as
(24.18) (implicit Euler method) In the compact notation mentioned above, we can also write this as
(24.19) (implicit Euler method) This is the same as the Euler method in Eq. (24.13), except that the derivative function is evaluated at the advanced time, f(y(t + ∆t), t + ∆t), instead of at the old time, f(y(t), t). The problem with this is that evolving y(t) to y(t + ∆t) requires knowing y(t + ∆t) already, which is of course why the method is implicit. There are a few ways to solve this equation at each time step. For example, the most straightforward is to use y(t) as a guess for y(t + ∆t), and plug it into f(y(t + ∆t), t + ∆t). Then generate the next guess for y(t + ∆t). Plug it back in, and keep guessing until the process converges to a solution, which is the one you want. This procedure is called fixed-point iteration, since the iteration converges to the desired steady state, or fixed point. Obviously, the implicit Euler method is a lot more work than the explicit counterpart, so why bother? What we gain in exchange for the extra effort is improved stability properties. Consider the simple model problem
(24.20) where \alpha is some large, positive constant. Obviously the solution is simply
(24.21) The (explicit) Euler method gives the update method
(24.22) Clearly, this method has rather serious problems if \alpha∆t > 1, since the solution y1 will be negative if y0 is positive, which definitely shouldn’t happen: the the analytic solution says the initially positive solution should stay that way. It’s pretty clear that in fact the solution will oscillate about zero, changing sign on each iteration. Things get even worse if \alpha∆t > 2, since now the coefficient of yn has a modulus of more than unity, and thus |yn| diverges to \infty. If \alpha is large, then it may take very small step sizes ∆t to obtain a stable recurrence. Even in more complicated systems of equations, the interesting dynamics may happen on relatively long time scales, but the step size may be limited to a short time interval by a fast decay (as in a set of rate equations with vastly different decay rates) to obtain a stable solution. Such a system of equations is called a stiff system. The implicit Euler method helps here, since the recursion now becomes
(24.23)
24.3.1 Second-Order Methods¶
24.3 Runge–Kutta Methods which when solved for yn+1 becomes
yn
(24.24) This recurrence is clearly stable (and not oscillatory) for any ∆t > 0, a distinct improvement over the explicit case. Of course, for a large step the solution still wouldn’t be very accurate. Often it’s the case that with an explicit method, a time step small enough to get an accurate solution also is small enough to guarantee stability. However, it sometimes helps to use an implicit method if a larger time step than would be explicitly stable gives an adequately accurate solution. The implicit Euler method is also called the backward Euler method, while the explicit Euler method is called the forward Euler method. Again, the improved stability properties inherited by the implicit method here also generally apply to more complicated implicit methods. Note that when going to the effort of using an implicit method, it is probably worth using a better method than the simple implicit Euler method. 24.3 Runge–Kutta Methods Now we will seek finite-difference methods for solving ODEs that are higher order, or that is to say methods that have global truncation errors that are O(∆tn) with n > 1. Of course, these methods will be more complicated, so again, why bother? Obviously, with a higher-order method, the step size ∆t needed to generate a solution with a particular accuracy will be larger than for a lower-order method. Often it is the case that the number of steps required for a high-order method is drastically smaller than for a low-order method, and there is correspondingly a large overall savings in computer time to obtain a particular level of accuracy. Furthermore, with a low-order method, it may not even be possible to achieve decent accuracy with a low-order method. In principle, by making the step size very small, you could achieve any desired accuracy with the Euler method. However, when the finite differences start becoming comparable to machine epsilon, machine rounding errors will dominate any truncation errors, providing an effective accuracy bound for a given numerical precision. A high-order method helps here because the step size to achieve a good accuracy will be much larger than for the low-order method. The accumulation of many small differences is replaced by a more complicated sequence of larger finite differences, which are less susceptible to roundoff errors. 24.3.1 Second-Order Methods To construct a second-order method, consider the following forward and backward Taylor expansions, ex- panded to higher order:
2! ¨y(t) + ∆t3 3! ...y (t) + ∆t4 4! y(4)(t) + O(∆t5)
2! ¨y(t) −∆t3 3! ...y (t) + ∆t4 4! y(4)(t) + O(∆t5) (24.25) Subtracting these, we find
3! ...y (t) + O(∆t5), (24.26) and now letting ∆t −\rightarrow t + ∆t and then ∆t −\rightarrow ∆t/2,
4 \cdot 3! ...y (t) + O(∆t5). (24.27) Thus, the finite-difference method
(24.28)
Chapter 24. Ordinary Differential Equations is O(∆t3) (i.e., locally second-order accurate). It also has what turns out to be the useful property that all even powers in the expansion of the error term vanish. This is because of the symmetry of the expression about the middle time tn+1/2. However, it requires that we know the intermediate solution yn+1/2, which of course we don’t know, and we wouldn’t even necessarily want to know it if we are trying to generate the solution samples yn. To deal with the intermediate value, we need to come up with an approximation for it. There are a number of ways to do this, each leading to different methods. One possibility is to take half of an Euler step to compute yn+1/2 and then use this result in the formula (24.28). Thus noting the two formulae
(24.29) we can write
(24.30) since the error in using d2 in place of yn+1/2 is O(∆t2), and with the extra factor of ∆t leads to an over- all O(∆t3) error. This leads to the midpoint method, an example of a second-order Runge–Kutta method:
2 f(yn, tn)
(24.31) (Runge–Kutta method, order 2) Again, the local truncation error is O(∆t3), while the global truncation error is O(∆t2). There are other Runge–Kutta methods that achieve second-order accuracy, which involve different choices for intermediate steps. 24.3.1.1 Variations The above choice of the second-order Runge–Kutta method is by no means unique.1 We can parameterize second-order Runge–Kutta methods more generally as
(general Runge–Kutta method, order 2) (24.32) That is dn is the result yn+a of taking an Euler step of a ∆t. The final solution yn+1 is a linear combination of taking a normal Euler step from (yn, tn) and an Euler step from (yn+a, tn+c). To see how we can obtain more second-order methods, we can write out the Taylor expansions of these evolution equations as
- b2 ∆t f(yn, tn) + b2 c ∆t2 \partial tf(yn, tn) + b2 a ∆t2 f(yn, tn) \partial yf(yn, tn) + O(∆t3). (24.33) Comparison to the exact Taylor expansion of yn+1
2 ¨yn + O(∆t3)
2 \partial tf(yn, tn) + ∆t2 2 f(yn, tn) \partial yf(yn, tn) + O(∆t3), (24.34) 1For a detailed treatment, see Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, (Cambridge, 1996).
24.3.2 Fourth-Order and General Schemes¶
24.3 Runge–Kutta Methods and matching terms up to O(∆t2) leads to the conditions
b2 c = 1 2, a = c. (order 2 Runge–Kutta conditions) (24.35) Any method of the form (24.32) with coefficients satisfying these three conditions is second order, which clearly leaves much wiggle room. Thus, we may regard the general order-two method as parameterized by the single parameter a:
1 −1 2a ∆t f(yn, tn) + 1 2a∆t f(dn, tn + a∆t). (general Runge–Kutta method, order 2) (24.36)
24.3.2 Fourth-Order and General Schemes To construct higher-order general schemes, recall that to evolve yn+1 from yn, we are performing the integral from tn to tn+1:
Z tn+1 tn dt′ f(y(t′), t′). (24.37) In constructing the second-order method above, we chose Z tn+1 tn
(24.38) as we showed in Eq. (24.27). The midpoint approximation is one of the simplest choices to represent the integral. We could of course sample the interval (tn, tn+1) with additional points, hoping to approximate the integral more exactly using a higher-resolution discrete sum, thus ending up with higher order methods. Thus, we can for example write
X \alpha
(24.39) where the coefficients b\alpha and c\alpha are to be determined to give the best possible approximation to the integral. Of course, more terms in the sum should lead to higher-order methods. And, of course, the problem with this scheme are that the solution at the intermediate times tn + c\alpha ∆t is not known, and thus must be estimated numerically. Following the idea from the second-order case, we write down approximations for each of the interme- diate solutions, and then compute a linear combination of the intermediate solutions to obtain an accurate
2Anthony Ralston, ‘‘Runge–Kutta Methods with Minimum Error Bounds,’’ Mathematics of Computation 16, 431 (1962).
Chapter 24. Ordinary Differential Equations write down the general form for the Runge–Kutta methods as d1 = yn
...
m−1 X \alpha=1
m X \alpha=1
(24.40) This set of formulae is commonly rewritten in the form
... km = ∆t f
yn + m−1 X \alpha=1 am\alpha, tn + cm ∆t !
m X \alpha=1
(general form for Runge–Kutta methods) (24.41) which is somewhat more efficient in terms of coding. The coefficients for fairly low-order systems can be chosen by matching the Taylor expansions of these formulae with the exact Taylor expansion, just as in the second-order case. This is cumbersome, and more elegant methods are available for constructing high-order methods.3 Here we will just quote some results at fourth order, where there is a good compromise between complexity of the method and accuracy. One very popular method is the ‘‘classical’’ fourth-order Runge–Kutta method, which takes four intermediate steps and can be written explicitly as
k2 = ∆t f yn + k1 2 , tn + ∆t k3 = ∆t f yn + k2 2 , tn + ∆t
6 + k2 3 + k3 3 + k4 6 . (classical fourth-order Runge–Kutta method) (24.42) Again, it is possible to choose the coefficients to not only cancel errors up to fourth order, but also to 3Arieh Iserles, op. cit.
24.3.3 Implicit Runge-Kutta Methods¶
24.3 Runge–Kutta Methods minimize the coefficient of the O(∆t5) error term.4¶
(fourth-order Runge–Kutta) (24.43) The coefficients for the time increments are the free parameters, and are given by c2 = 2 \approx 0.4000 0000 0000 0000 c3 = 7 8 −3 \sqrt \approx 0.4557 3725 4218 7894, (fourth-order Runge–Kutta coefficients) (24.44) while the coefficients for the linear combination are
\sqrt \approx 0.1747 6028 2262 6904
\sqrt 5) \approx −0.5514 8066 2878 7329 b3 = \sqrt 5 −10038 \approx 1.2055 3559 9396 5235 b4 =
\sqrt \approx 0.1711 8478 1219 5190, (fourth-order Runge–Kutta coefficients) (24.45) and the coefficients for intermediate linear combinations of solutions are a21 = 2 \approx 0.4000 0000 0000 0000
\sqrt 5 −963) \approx 0.2969 7760 9247 7536
\sqrt 5) \approx 0.1587 5964 4971 0358 a41 = 2094 \sqrt 5 −3365 \approx 0.2181 0038 8225 9205
\sqrt \approx −3.0509 6514 8692 9308
\sqrt 5) 240845 \approx 3.8328 6476 0467 0103. (fourth-order Runge–Kutta coefficients) (24.46) This method produces an error bound about half that of the classical formula (24.42), though it takes somewhat more effort to program and a few extra arithmetic operations. 24.3.3 Implicit Runge–Kutta Methods The Runge–Kutta methods we have considered so far are explicit. Looking for example at Eq. (24.40), the mth intermediate solution dm only depends on the previous dj where j < m. Of course this constraint can 4Anthony Ralson, op. cit.
24.3.4 Adaptive Stepping¶
Chapter 24. Ordinary Differential Equations be relaxed, and each dm can depend on function evaluations at every other dm. This generalization leads to implicit Runge–Kutta (Gauss–Legendre) methods. For example, the two-stage scheme
(two-stage implicit Runge–Kutta method) (24.47) with coefficients a11 = 1 a12 = 1 4 − \sqrt a21 = 1 4 + \sqrt a22 = 1 b1 = 1 b2 = 1 c1 = 1 2 − \sqrt c2 = 1 2 + \sqrt (fourth-order implicit Runge–Kutta coefficients) (24.48) turns out to have fourth-order accuracy.5 This seems more compact than the four-stage explicit Runge–Kutta methods above, but of course there is extra complexity associated with functional iteration or whatever is used to deal with the implicit equations. However, the implicit method should be more stable and robust than the explicit counterpart, and better able to deal with stiff problems. 24.3.4 Adaptive Stepping An important concept in numerical integration is adaptive time stepping. The idea is that, along with estimating the solution, we can also estimate the error in the solution. Then, assuming we have some goal to achieve for the local accuracy (truncation error), or tolerance, we can perform a sort of feedback to the step size to just meet the specified tolerance. In other words, if the error is too large, then the step is rejected, and a smaller step is taken. If the error estimate is below tolerance, then the next time step is chosen to be larger, so that the desired accuracy is achieved with the minimum level of work. This strategy has two main advantages. First, it removes the problem of choosing a sensible step size from the user, since the algorithm handles this automatically. This reduces potential for user error and wasted cpu cycles. More importantly, it may be that in some solutions, the nature of the solution may be continually in flux, with long steps appropriate during slowly varying intervals, and short steps necessary during periods of intense action. If the quiet intervals are large, the time savings gained by using an adaptive method can be substantial. In the case of Runge–Kutta methods, the methods we have outlined don’t already have an error estimate built in. Thus, we must do extra work. For example, after computing yn+1 from yn by one Runge– Kutta step, we could recompute it using two half-steps, using the same algorithm. This is the strategy used, for example, in Richardson extrapolation (Section 26.3). However, the more common method used with Runge–Kutta integrators is to use two different methods of different order. A very popular choice is called the (explicit) Runge–Kutta–Fehlberg method,6 where a fourth-order and a fifth-order method are computed simultaneously on each time step, and the difference between the two acts as an estimate for the O(∆t5) truncation error of the fourth-order method. The Runge–Kutta–Fehlberg method cleverly arranges both methods to make use of the same set of function evaluations, to minimize the extra effort beyond the basic fourth-order method. Such methods are called embedded, because for example the fourth- order method is ‘‘embedded’’ within the fifth-order method. Because we know how the error scales, it is reasonably straightforward to decide how to rescale the step size to achieve the desired error tolerance. 5Arieh Iserles, op. cit. 6Erwin Fehlberg, ‘‘Low-Order Classical Runge-Kutta Formulas with Stepsize Control and their Application to some Heat Transfer Problems,’’ NASA Technical Report R-315 (1969), Table III. See also Arieh Iserles, op. cit., p. 84.
24.4 Multistep and Predictor–Corrector Methods Adjusting the time step is fairly straightforward with explicit Runge–Kutta methods, because the solution yn+1 only depends on yn. This property does not hold for the predictor-corrector methods below, and they are concomitantly much more complicated to code as adaptive-step methods. The fixed-step Runge–Kutta methods are fairly straightforward to code, but the extra overhead in- volved with the adaptive methods makes them slightly more complicated to code. Fortunately, there are free, high-quality implementations available, such as RKSUITE on Netlib.7 Now the adaptive-step method will obviously generate a solution that is not uniformly sampled in time, whereas often a uniform sampling is desirable. There are two main ways to deal with this. One is to go ahead and generate the solution with nonuniform samples, and then use numerical interpolation (e.g., cubic splines) to transfer the solution to the desired set of sample times. Another popular method is to decide on the desired sample times in advance, and then call the integrator to evolve the solution forward from one sample time to the next. The integrator steps adaptively to the next time, and terminates at the desired output time by adjusting its final step to exactly hit the output time. The integrator subroutine is thus called once for each output time step (whereas in the former method, the integrator subroutine is called only once for the whole solution). The latter method is especially common in compiled languages (e.g., Fortran or C), as opposed to ‘‘canned’’ languages like Octave. 24.4 Multistep and Predictor–Corrector Methods Another important class of higher-order numerical methods goes by the name of multistep methods. The idea is that, as in the Runge–Kutta methods, to achieve higher order you need to estimate the solution at multiple points in time. Runge–Kutta methods do this by estimating the solution at several intermediate times between t and t+∆t. Multistep methods, however, use the additional points in the past, say at t−∆t, t −2∆t, and so on, to achieve the higher order. We won’t go into the derivation of these methods here,8 but for example, the second-order Adams–Bashforth method is
3 2f(yn, tn) −1 2f(yn−1, tn−2) , (24.49) while the third-order method is
23 12f(yn, tn) −4 3f(yn−1, tn−2) + 5 12f(yn−2, tn−2) . (24.50) Common implicit versions of these methods are Adams–Moulton methods and backward differentia- tion formulae, which are better when improved stability is needed or when solving stiff problems. Implementations of these methods can be complicated, especially in implementing adaptive-step ver- sions. Since the next step may rely on several steps in the past, changing future time steps require either interpolating the past solution onto the new time grid or using much more complicated formulae that can handle nonuniform time steps. Further, imagine: how would you start a high-order multistep method, given only initial data y(t0)? You could, for example step backwards using a low-order method (e.g., Euler) to generate the past steps. Sophisticated implementations start with the Euler method, and increase the or- der as the solution proceeds, adaptively choosing both the order of the method and the stepsize needed to maximize efficiency at the desired error tolerance. Multistep methods find their biggest utility in predictor–corrector methods. The basic idea is to use two methods, one explicit and one implicit, such as an Adams–Bashforth and an Adams–Moulton method of the same order, and run them together. The explicit method (the predictor), generates the updated solution at the future time, which is then used as a starting point for iterating the implicit method (the corrector) to convergence. The difference between the two, or perhaps whether or not the iteration converges quickly enough, is used to determine whether or not the current step size is adequate. There are a number of different 7http://www.netlib.org/ode/rksuite/ 8See Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, (Cambridge, 1996), Chapter 2.
Chapter 24. Ordinary Differential Equations strategies here, but these methods are again quite complicated to implement. Fortunately, high-quality, open- source integrators are available, such as the venerable ODEPACK9 and others10 at Netlib. These methods have the advantage of efficiently producing solutions of very high accuracy. However, they perform best on very smooth solutions. In cases where certain errors are involved in computing the derivative function (as when lookup tables are used), these methods may not perform so well. They may also not perform well in cases where the integrator must be restarted often, as when the solution has discontinuities (as in quantum-jump trajectories), since the restarting overhead can be substantial. In these cases it better to fall back to the less elegant but more robust Runge–Kutta methods. 24.5 Exercises Problem 24.1 An alternate definition of the limit of a function is as follows. We say y = lim x\rightarrow x0 f(x) (24.51) if for every ϵ > 0, there is a \delta > 0 such that if |x −x0| < \delta then |f(x) −f(xn)| < ϵ. Essentially, this is saying that a sufficiently small perturbation to x can make an arbitrarily small perturbation to f(x). Show that this ‘‘ϵ-\delta’’ definition of the limit is equivalent to the definition in terms of limits of sequences. 9http://www.netlib.org/odepack/ 10http://www.netlib.org/ode/