Euler Method – How to Create Adaptive Step Size

Question

I think Euler's Method is a great way to simulate ODE:s. It's not the most accurate, but it's the fastest and simplest.

Euler's Method is usually used with fixed step size, where $k$ is the step size larger than $0$ and $\dot x = f(x,u)$ is our ODE function.

To simulate forward Euler, just iterate this equation:

$$x_{i+1} = x_i + k f(x_i, u)$$

To improve stability for Euler's method, then the step size $k$ needs to be adaptive. The smaller choice of the step $k$, the more stability is guaranteed, but the less accurate the simulation will become. Instead, the adaptive method for the step size is often used to find the best $k$.

Question:

I can I improve stability for Euler's method by implementing adaptive step size?

Answer 1

You will need a second order method to compute the error estimate. You can use any of \begin{array}{c|cc} 0&\\ c&c\\ \hline &1-\frac1{2c}&\frac1{2c} \end{array} The Heun method, $c=1$, has the advantage that one can re-use the second slope as the first slope in the next step.

h=h0
k1 = f(t,y)    # derivative to be used for the explicit Euler step
while t<tf:
    k2 = f(t+h,y+h*k1)
    err = 0.5*abs(k1-k2)
    fac = tol/err  # better use relative error
    if fac > 1: # accept step
        y += h*k1; k1 = k2  # Euler step, derivative was already computed
        t += h
        # collect t,y in some output structure or print to file
    h = 0.9*fac*h
    # for exploration collect t,h in output structure or a second file
# return output structure

In the comments the question for a variable-step Heun method came up. I cited the third-order methods in Butcher: "Low-order methods", slide 26, as any explicit 3rd order method contains second order method in its first two stages. In Feldman: "Variable step size methods" you can find several tested low-order embedded methods. For instance one based on Heun's method attributed to Fehlberg that probably started him to search for higher order embedded methods \begin{array}{c|ccc} 0&\\ 1&1\\ \frac12&\frac14&\frac14\\ \hline &\frac16&\frac16&\frac46 \end{array} This computes from the Heun data of the first two stages the point and derivative at the midpoint of the interval. Then applying the Simpson rule to these approximations results in an order 3 approximation.

h=h0
while t<tf:
    k1 = f(t,y)
    k2 = f(t+h,y+h*k1)
    k12 = 0.5*(k1+k2)   # combined slope for Heun step
    k3 = f(t+0.5*h,y+0.5*k12)
    err = (2/3)*abs(k12-k3)
    fac = tol/err  # better use relative error
    if fac > 1: # accept step
        y += h*k12  # Heun step, this is a new point
        t += h
        # collect t,y in some output structure or print to file
    h = 0.9*fac*h
    # for exploration collect t,h in output structure or a second file
# return output structure