Improved Euler and Second-Order Runge-Kutta, differential equations lesson

Euler’s method reads the slope only at the start of a step, so it cannot see how the slope changes before the step ends. If one more slope is read somewhere inside the step, where should it be read, and how should the two readings be combined? Throughout, $t$ is the independent variable, $y' = dy/dt$ , the problem is $y' = f(t, y)$ with $y(t_0) = y_0$ , and the grid is $t_n = t_0 + nh$ with step size $h > 0$ . This lesson works through two answers, Heun’s method and the midpoint method, explains why both have global error of order $h^2$ , and shows that they are still different methods.

Key terms

Stage slopes

Definition An explicit Runge-Kutta step from $(t_n, y_n)$ with $s$ stages evaluates $f$ at $s$ trial points across the step. Stage $i$ has the time $t_n + c_i h$ and a stage state $\tilde y_i$ built from $y_n$ and the earlier stage slopes, and its stage slope $k_i$ is the value of $f$ there,

\begin{aligned} \tilde y_i &= y_n + h \sum_{j=1}^{i-1} a_{ij} k_j, \\[0.4em] k_i &= f(t_n + c_i h, \tilde y_i), \\[0.4em] y_{n+1} &= y_n + h \sum_{i=1}^{s} b_i k_i. \end{aligned}

The first stage has an empty sum, so $\tilde y_1 = y_n$ , and the methods in this unit also take $c_1 = 0$ , so the first slope is read at the current point. Each $k_i$ is a slope, in units of $y$ per unit of $t$ , so every change in $y$ , whether in a stage state or in the update, is $h$ times a combination of slopes. The stage states are trial values used only to evaluate slopes. They are not points of the numerical solution, and only $y_{n+1}$ is kept. Each stage costs one evaluation of $f$ .

In plain words Each stage is a trial point inside the step where the method reads the slope field. The method then moves from $y_n$ by a weighted combination of those slopes and throws the trial points away.

Example Take Heun’s method, described below, on $y' = y$ with $y(0) = 1$ and $h = \tfrac15$ . Stage one is read at $(0, 1)$ and gives $k_1 = 1$ . Stage two is read at the predicted point $(\tfrac15, \tfrac65)$ and gives $k_2 = \tfrac65$ . The accepted value at $t = \tfrac15$ is $\tfrac{61}{50} = 1.22$ , not the stage state $\tfrac65 = 1.2$ . The midpoint method reads its second stage at $(\tfrac1{10}, \tfrac{11}{10})$ instead.

Heun’s method

Definition Heun’s method is the two-stage method

\begin{aligned} k_1 &= f(t_n, y_n), \\[0.4em] k_2 &= f(t_n + h, y_n + h k_1), \\[0.4em] y_{n+1} &= y_n + \frac{h}{2}(k_1 + k_2). \end{aligned}

The second stage is read at the end of the step, at the state predicted by an Euler step, and the update uses the average of the two slopes. The Euler prediction $y_n + hk_1$ only locates $k_2$ and is not the result of the step. Many textbooks call this method the improved Euler method, while some attach the names improved or modified Euler to other formulas, including the midpoint method, so a name should always come with its formula. When $f$ depends on $t$ alone, the step is the trapezoidal rule for the integral of $f$ over the step.

In plain words Take a trial Euler step, read the slope where it lands, and then redo the step from the start with the average of the starting slope and that endpoint slope.

Example Take $y' = y$ , $y(0) = 1$ and $h = \tfrac15$ . In the first step $k_1 = 1$ , the prediction is $\tfrac65$ , so $k_2 = \tfrac65$ , the average slope is $\tfrac{11}{10}$ , and

y_1 = 1 + \tfrac15 \cdot \tfrac{11}{10} = \tfrac{61}{50}.

In the second step, from $t_1 = \tfrac15$ , the slopes are $k_1 = \tfrac{61}{50}$ and $k_2 = \tfrac{61}{50} \cdot \tfrac65 = \tfrac{183}{125}$ . Their average is $\tfrac{671}{500}$ , so

\begin{aligned} y_2 &= \tfrac{61}{50} + \tfrac15 \cdot \tfrac{671}{500} \\[0.4em] &= \tfrac{3721}{2500} = 1.4884. \end{aligned}

The exact solution $e^t$ , which equals its own derivative and starts at $1$ , gives $e^{0.4} \approx 1.49182$ , so the error is about $-0.00342$ . Euler’s method with the same step gives $1.2^2 = 1.44$ , with error about $-0.0518$ .

Midpoint method

Definition The explicit midpoint method is the two-stage method

\begin{aligned} k_1 &= f(t_n, y_n), \\[0.4em] k_2 &= f\bigl(t_n + \tfrac{h}{2}, y_n + \tfrac{h}{2} k_1\bigr), \\[0.4em] y_{n+1} &= y_n + h k_2. \end{aligned}

The first slope only predicts the state at the middle of the step, and the single slope read there is used for the whole step. When $f$ depends on $t$ alone, the step is the midpoint rule for the integral of $f$ over the step. The midpoint method and Heun’s method have the same order but different stages, and in general they give different numbers.

In plain words Take half an Euler step to estimate where the solution is at the middle of the step, read the slope there, and use that one slope from the start of the step to its end.

Example On $y' = y$ with $y(0) = 1$ and $h = \tfrac15$ , the midpoint state is $1 + \tfrac1{10} = \tfrac{11}{10}$ , so $k_2 = \tfrac{11}{10}$ and $y_1 = 1 + \tfrac15 \cdot \tfrac{11}{10} = \tfrac{61}{50}$ , the same value as Heun’s method. That agreement does not hold for every equation. On $y' = t^2$ with $y(0) = 0$ , one step of size $h$ gives

\begin{gathered} \text{Heun}\quad \tfrac{h}{2}\bigl(0 + h^2\bigr) = \tfrac{h^3}{2}, \\[0.4em] \text{midpoint}\quad h\bigl(\tfrac{h}{2}\bigr)^2 = \tfrac{h^3}{4}, \end{gathered}

while the exact value is $y(h) = h^3/3$ . The one-step errors, computed minus exact, are $h^3/6$ and $-h^3/12$ . Both are of order $h^3$ , but the two methods sample the slope at different times and give different answers. With $h = \tfrac15$ the three numbers are $\tfrac1{250}$ , $\tfrac1{500}$ and $\tfrac1{375}$ .

Second-order Runge-Kutta methods

Definition Consider the two-stage explicit methods with the stage slopes

\begin{aligned} k_1 &= f(t_n, y_n), \\[0.4em] k_2 &= f(t_n + ch, y_n + ch\,k_1) \end{aligned}

and the update

y_{n+1} = y_n + h(b_1 k_1 + b_2 k_2).

Along a solution the chain rule gives $y'' = f_t + f_y f$ , so for smooth $f$ , with $f$ and its partial derivatives evaluated at $(t_n, y(t_n))$ ,

\begin{aligned} y(t_{n+1}) &= y(t_n) + h f \\[0.4em] &\quad + \tfrac{h^2}{2}(f_t + f_y f) \\[0.4em] &\quad + O(h^3). \end{aligned}

A step started from $y(t_n)$ has $k_2 = f + ch(f_t + f_y f) + O(h^2)$ , so it reproduces these terms for every smooth $f$ exactly when

\begin{gathered} b_1 + b_2 = 1, \\[0.4em] b_2 c = \tfrac12. \end{gathered}

The methods that satisfy both conditions are the second-order Runge-Kutta methods of this form. The second condition forces $c \ne 0$ and leaves one free parameter, with $b_2 = 1/(2c)$ and $b_1 = 1 - 1/(2c)$ . Heun’s method has $c = 1$ and $b_1 = b_2 = \tfrac12$ , the midpoint method has $c = \tfrac12$ , $b_1 = 0$ and $b_2 = 1$ , and Ralston’s method has $c = \tfrac23$ , $b_1 = \tfrac14$ and $b_2 = \tfrac34$ . No two-stage explicit method has order $3$ for every smooth $f$ .

In plain words The second slope carries information about how the slope changes during the step, which is the $h^2$ term of the Taylor series that Euler’s method misses. The two conditions say that the weights add up to one and that the second stage is placed and weighted so that this $h^2$ term comes out right. Any nonzero placement $c$ works once its weight is adjusted to match.

Example Heun’s method satisfies $\tfrac12 + \tfrac12 = 1$ and $\tfrac12 \cdot 1 = \tfrac12$ , and the midpoint method satisfies $0 + 1 = 1$ and $1 \cdot \tfrac12 = \tfrac12$ . On $y' = y$ every method of the family multiplies the value in each step by

R(h) = 1 + h + \tfrac{h^2}{2},

which matches $e^h$ through $h^2$ , and $e^h - R(h) = h^3/6 + O(h^4)$ . With $h = 0.2$ the difference is $0.001403$ , close to $h^3/6 \approx 0.001333$ . Compare a formula that reads $k_2$ at the end of the step, as Heun’s method does, but advances with $y_{n+1} = y_n + hk_2$ . It has $b_2 c = 1 \ne \tfrac12$ , its factor on $y' = y$ is $1 + h + h^2$ , and its one-step error is of order $h^2$ , so it is only first order.

Second-order global error

Definition Let $f$ have continuous partial derivatives through order two on an open set containing the graph of the exact solution over $[t_0, T]$ , a closed interval inside its interval of existence. Then for every second-order method of the family above there are constants $C$ and $h_0 > 0$ such that

\lvert y_n - y(t_n) \rvert \le C h^2

for every grid time $t_n \le T$ whenever $0 < h \le h_0$ . The local truncation error, the error of one step started from the exact value, is of order $h^3$ , and with the normalized convention that divides by $h$ it is of order $h^2$ . In a typical problem, once $h$ is small, halving the step divides the global error by about $4$ . The smoothness matters. If $f$ jumps in $t$ at a time inside a step, that one step can have an error of order $h$ , and the global error can then be only of order $h$ . When $f$ is smooth on each side of the jump, stopping the computation at the jump and restarting there with the new formula restores second order. The order also describes small steps only. On $y' = -ay$ with $a > 0$ , every method of the family multiplies by $1 - ah + (ah)^2/2$ per step, and the computed values decay exactly when $0 < h < 2/a$ .

In plain words With two slopes per step, the error falls like $h^2$ . Halving the step doubles the work and, for a smooth problem with a small enough step, makes the error about four times smaller. A second-order label promises nothing about a step that is too large.

Example Take $y' = y^2$ with $y(0) = 1$ , whose solution $y = 1/(1 - t)$ exists on $(-\infty, 1)$ , and compute up to $t = \tfrac12$ , where $y = 2$ . The signed errors at $t = \tfrac12$ are

$h$	Heun	Midpoint
$0.1$	$-1.670 \times 10^{-2}$	$-2.294 \times 10^{-2}$
$0.05$	$-4.598 \times 10^{-3}$	$-6.579 \times 10^{-3}$
$0.025$	$-1.201 \times 10^{-3}$	$-1.758 \times 10^{-3}$
$0.0125$	$-3.065 \times 10^{-4}$	$-4.541 \times 10^{-4}$

Each halving divides the errors by a factor that approaches $4$ , from $3.6$ to $3.9$ for Heun’s method and from $3.5$ to $3.9$ for the midpoint method. Both methods are second order, yet here the midpoint errors are about $1.4$ to $1.5$ times larger, so the same order does not mean the same accuracy.

Check yourself

Take one step of Heun’s method and one step of the midpoint method, each with $h = 0.2$ , for $y' = -y^2$ , $y(0) = 1$ . Compare both with the exact solution $y = 1/(1 + t)$ .

Show the answer

Both start with $k_1 = f(0, 1) = -1$ . Heun’s method predicts $1 + 0.2(-1) = 0.8$ at $t = 0.2$ , reads $k_2 = -(0.8)^2 = -0.64$ there, and gives $y_1 = 1 + 0.1(-1 - 0.64) = 0.836$ . The midpoint method predicts $1 + 0.1(-1) = 0.9$ at $t = 0.1$ , reads $k_2 = -0.81$ , and gives $y_1 = 1 + 0.2(-0.81) = 0.838$ . The proposed exact solution satisfies $y' = -1/(1 + t)^2 = -y^2$ and $y(0) = 1$ , and $y(0.2) = 1/1.2 \approx 0.83333$ . The signed errors are about $0.00267$ for Heun’s method and $0.00467$ for the midpoint method, compared with about $-0.0333$ for one Euler step, which gives $0.8$ .

A classmate argues that Heun’s method is second order, so the step $h = 0.25$ must be fine for $y' = -10y$ , $y(0) = 1$ . What does Heun’s method actually do here, and which steps make its values decay?

Show the answer

Each step multiplies by $1 + z + z^2/2$ with $z = -10h$ , that is, by $1 - 10h + 50h^2$ . The step $h = 0.25$ gives the factor $1 - 2.5 + 3.125 = 1.625$ , so the values $1$ , $1.625$ , $2.641$ and so on grow, while the exact solution $e^{-10t}$ decays. With $x = 10h$ the factor is $1 - x + x^2/2$ , which is never below $\tfrac12$ , so decay requires only that it be below $1$ . That means $x^2/2 < x$ , or $0 < h < 0.2$ . With $h = 0.1$ the factor is $0.5$ and the values decay. The order describes how the error shrinks as $h \to 0$ , and it says nothing about whether one particular step is too large.

For $y' = t - y$ , $y(0) = 1$ and $h = 0.2$ , Heun’s method and the midpoint method both give $y_1 = 0.84$ . Check this, and explain why the two methods agree on this equation but not on $y' = -y^2$ from question 1.

Show the answer

Heun’s method has $k_1 = -1$ , the prediction $0.8$ at $t = 0.2$ and $k_2 = 0.2 - 0.8 = -0.6$ , so $y_1 = 1 + 0.1(-1.6) = 0.84$ . The midpoint method has the state $0.9$ at $t = 0.1$ and $k_2 = 0.1 - 0.9 = -0.8$ , so $y_1 = 1 + 0.2(-0.8) = 0.84$ . In general, if $f(t, y) = \alpha t + \beta y + \gamma$ with constants $\alpha$ , $\beta$ and $\gamma$ , then $k_2 = k_1 + ch(\alpha + \beta k_1)$ exactly, and since $b_1 + b_2 = 1$ the update is

\begin{aligned} y_{n+1} &= y_n + h k_1 \\[0.4em] &\quad + b_2 c\,h^2(\alpha + \beta k_1). \end{aligned}

Every second-order method has $b_2 c = \tfrac12$ , so all of them take exactly the same step on such an equation. For a nonlinear right side such as $-y^2$ , the second slope contains further terms in $ch$ that depend on $c$ , and the methods differ. For comparison, the exact solution $y = t - 1 + 2e^{-t}$ gives $y(0.2) \approx 0.83746$ .

Improved Euler and Second-Order Runge-Kutta

Key terms

Stage slopes

Heun’s method

Midpoint method

Second-order Runge-Kutta methods

Second-order global error

Check yourself

How this lesson connects

Builds on

Used later in