Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I was recently watching a tutorial on Euler's method for approximating differential equations, and the whole time I was thinking "why can't you just take the limit of the step size $h$ as it goes to $0$, and get an exact or near-exact approximation of the differential equation?" So that is my question: is it possible to do that?

share|improve this question
    
In theory, yes. In practice, rounding error accumulates and this doesn't work. –  Michael Lugo Oct 28 '11 at 15:22
2  
"exact approximation" is oxymoronic, yes? –  J. M. Oct 28 '11 at 15:28
1  
It is a nice idea. For an explicit solution of the DE, we would need to compute the limit explicitly. Even for nice functions, this can be quite hard. Works well for a few cases, for example $\frac{dy}{dx}=ky$, but that's pretty easy anyway. The same issue arises with the Riemann integral. It is defined as the limit of a certain sum, but apart from a few introductory examples, that's not how we find exact expressions for definite integrals, because finding the limit is too hard. Instead, we prove a result that connects the Riemann integral to antiderivatives. –  André Nicolas Oct 28 '11 at 15:38

1 Answer 1

up vote 6 down vote accepted

Not really. Recall that the Euler method for $y^\prime=f(x,y)$ takes the form

$$\begin{align*}x_{k+1}&=x_k+h\\y_{k+1}&=y_k+hf(x_k,y_k)\end{align*}$$

for some stepsize $h$. Taking $h=0$ here is equivalent to not moving at all!

However, your idea can be made slightly more practical. Consider the application of the Euler method with stepsize $h/2$:

$$\begin{align*}x_{k+1/2}&=x_k+h/2\\y_{k+1/2}&=y_k+hf(x_k,y_k)/2\\x_{k+1}&=x_{k+1/2}+h/2\\y_{k+1}&=y_{k+1/2}+hf(x_{k+1/2},y_{k+1/2})/2\end{align*}$$

The value of $y_{k+1}$ from these two steps can usually be expected to be a bit more accurate that the value of $y_{k+1}$ from the $h$-step Euler method. We can keep playing the game, taking 4 steps with stepsize $h/4$, 8 steps with stepsize $h/8$, and so on, yielding a sequence of estimates for $y_{k+1}$ corresponding to decreasing $h$.

One way one might estimate the result of what happens when $h\to 0$ is to take all those estimates of $y_{k+1}$ along with the associated stepsizes, and then fit an interpolating polynomial to those. For example, taking $y_{k+1}^{(0)}$ to be the result for stepsize $h$, $y_{k+1}^{(1)}$ the corresponding result for stepsize $h/2$, and $y_{k+1}^{(2)}$ the result for stepsize $h/4$, one can fit a quadratic interpolating polynomial to the three points $\{(h,y_{k+1}^{(0)}),(h/2,y_{k+1}^{(1)}),(h/4,y_{k+1}^{(2)})\}$, and then estimate the limit as $h\to 0$ by evaluating the interpolating polynomial thus obtained at 0.

"Lady Windermere's fan"

This scheme of using the interpolating polynomial to estimate the limit to 0 is called Richardson extrapolation. In practice, one certainly takes more than three points for the interpolation, and the order of the interpolating polynomial needed is estimated based on the behavior at past points. The idea of using Richardson extrapolation is due to Roland Bulirsch and Josef Stoer. The Bulirsch-Stoer method discussed here uses a slightly different method for integrating the differential equation (the modified midpoint method) from which the extrapolations are built (as well as a slightly modified extrapolation method), but is essentially the same idea as presented here.


Here's a tiny Mathematica demonstration of the Bulirsch-Stoer idea:

Table[(Apply[List, 
     y /. First[
       NDSolve[{y'[x] == x - y[x], y[0] == 1}, y, {x, 0, 1}, 
        MaxStepFraction -> 1, Method -> "ExplicitEuler", 
        InterpolationOrder -> 1, StartingStepSize -> 2^-k]]])[[4, 
    3, -2]], {k, 4}] - 2/E

{-0.235759, -0.102946, -0.0485411, -0.0236106}

Here we tried the Euler method on the differential equation $y^\prime=x-y$ with initial condition $y(0)=1$ for the stepsizes $h=1/2,1/4,1/8,1/16$, and compared the result at $x=1$ with the true value $y(1)=2/e$. As you can see the accuracy isn't too good for any of these.

Here's what happens after Richardson extrapolation using those same results from Euler:

InterpolatingPolynomial[
  Table[{2^-k, (Apply[List, 
       y /. First[
         NDSolve[{y'[x] == x - y[x], y[0] == 1}, y, {x, 0, 1}, 
          MaxStepFraction -> 1, Method -> "ExplicitEuler", 
          InterpolationOrder -> 1, StartingStepSize -> 2^-k]]])[[4, 
      3, -2]]}, {k, 4}], 0] - 2/E

0.0000823142

We end up with a result good to three or so digits, much better than even the result of Euler corresponding to $h=1/256$. Pretty good, I would say, considering that only $2+4+8+16=30$ evaluations of $f(x,y)=x-y$ were needed for a result with three-digit accuracy, while Euler with $256$ steps (and thus $256$ evaluations of $f(x,y)$) can only manage two good digits.

As an aside, Mathematica implements Bulirsch-Stoer internally in NDSolve[], as the option Method -> "Extrapolation".

share|improve this answer
    
thank you so much for this extremely in depth explanation, I really appreciate you time, and you answer. it really clarifies things for me. –  AlexW.H.B. Oct 30 '11 at 1:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.