1 Fundamental Methods of Extrapolation Fundamental Methods of Numerical Extrapolation With Applications Eric Hung-Lin Liu Keywords: numerical analysis, extrapolation, richardson, romberg, numerical differentiation, numerical integration Abstract Extrapolation is an incredibly powerful technique for increasing speed and ac curacy in various numerical tasks in scientific computing. As we will see, extrapolation can transform even the most mundane of algorithms (such as the Trap e zoid Rule) into an extremely fast and accurate algorithm, increasing the rate of conver gence by more than one order of magnitude. Within this paper, we will first present some background theory to motivate the derivation of Richardson and Romberg Extrapolation and provide support for their validity and stability as numerical techniques. Then we will present examples from differentiation and integration to demonstrate the incredible power of extrapolation. We will also give some MATLAB code to show some simple implentation methods in addition to discussing error bounds and differentiating between ”‘true”’ and ”‘false”’ convergence. 1 Introduction Extrapolation is an extremely powerful tool available to numerical analysts for improving the performance of a wide variety of mathematical methods. Since many “interesting” problems cannot be solved analytically, we often have to use computers to derive approximate solutions. In order to make this process efficient, we have to think carefully about our implementations: this is where extrapolation comes in. We will examine some examples in greater detail later, but we will briefly survey them now. First, consider numeric differentiation: in theory, it is performed by tending the step-size to zero. The problem is that modern computers have limited precision; as we will see, decreasing the step-size p ast 10−8 severely worsens the performance because of floating point round-off errors. The choice of step-size in integration is also important but for a different reason. Having a tiny step-size means more accurate results, but it also requires more function evaluations–this could be very expensive depending on the application at hand. The large numbers of floating point operations associated with a small step-size h as the negative effect of magnifying the problems of limited precision ( but it is not nearly as serious as with differentiation). Another example is the evaluation of (infinite) series. For slowly converging series, we will need to sum a very large number of terms to get an accurate result; the problem is that this could be very slow and again limited precision comes into play. Extrapolation attempts to rectify these problems by computing “weighted averages” ofX 2 Liu relatively poor estimates to eliminate error modes–specifically, we will be examining various forms of Richardson Extrapolation. Richardson and Richardson-like schemes are but a subset of all extrapolation methods, but they are among the most common. They have numerous applications and are not terribly difficult to understand and implement. We will demonstrate that various numerical methods taught in introductory calculus classes are in fact insufficient, but they can all be greatly improved through proper application of extrapolation techniques. 1.1 Background: A Few Words on Series Expansions We are very often able to associate some infinite sequence {An} with a known function A(h). For example, A(h) may represent a first divided differences scheme for differentiation, where h is the step-size. Note that h may be continuous or discrete. The sequence {An} is described by: An = A(hn), n ∈ N for a monotonically decreasing sequence {hn} that converges to 0 in the limit. Thus limh→0+ A(h) = A and similarly, limn→∞ An = A. Computing A(0) is exactly what we want to do–continuing the differentiation example, roughly speaking A(0) is the ”point” where the divided differences ”becomes” the derivative. Interestingly enough, although we require that A(y) have an asymptotic series expansion, we need not determine it uniquely. In fact we will see that knowing the pn is sufficient. Specifically, we will examine functions with expansions that take this form: s A(h) = A + αkhpk + O(hps+1 ) as h 0+ (1) →k=1 where pn 6= 0 ∀n, p1 < p2 < ... < ps+1, and all the αk are independent of y. The pn can in fact be imaginary, in which case we compare only the real parts. Clearly, having p1 > 0 guarantees the existence of limh→0+ A(h) = A. In the case where that limit does not exist, A is the antilimit of A(h), and then we know that at pi ≤ 0 for at least i = 1. However, as long as the sequence {pn} is monotonically increasing such that limk→∞ pk = +∞, (1) exists and A(h) P∞has the expansion A(h) ∼ A + k=1 αkhpk . We do not need to declare any conditions on the convergence of the infinite series. We do assume that the hk are known, but we do not need to know the αk, since as we will see, these constant factors are not significant. Now suppose that p1 > 0 (hence it is true ∀p), then when h is sufficiently small, A(y) approximates A with error A(h) − A = O(hp1 ). In most cases, p1 is not particularly small (i.e. p1 < 4), so it would appear that we will need to decrease h to get useful results. But as previously noted, this is almost always prohibitive in terms of round-off errors and/or function evaluations. We should not despair, because we are now standing on the doorstep of the fundamental idea behind Richardson Extrapolation. That idea is to use equation (1) to somehow eliminate the hp1 error term, thus obtaining a new approximation A1(h) −A = O(hh2 ). It is obvious that | A1(h) −A |<| A(h) −A | since p2 > p1. As we will see later, this is accomplished