The 4th-order Runge Kutta method for solving IVPs is to Heun's method as Simpson's rule is to the trapezoidal rule. It samples the slope at intermediate points as well as the end points to find a good average of the slope across the interval.

# Background

Useful background for this topic includes:

- 7. Taylor Series
- 14.2 Heun's Method

# References

- Bradie, Section 7.4, Runge-Kutta Methods, p.574.
- Mathews, Section 9.5, Runge-Kutta Methods, p.489.
- Weisstein, http://mathworld.wolfram.com/Runge-KuttaMethod.html.

# Interactive Maplet

# Theory

By using a similar strategy to the trapezoidal rule to find a better approximation to an IVP in Heun's method, consider now Simpson's rule, where not only the end points, but also the interior points of the interval are sampled.

The 4th-order Runge Kutta method is similar to Simpson's rule: a sample of the slope is made at the mid-point on the interval, as well as the end points, and a weighted average is taken, placing more weight on the slope at the mid point.

It should be noted that Runge-Kutta refers to an entire class of IVP solvers, which includes Euler's method and Heun's method. We are looking at one particularly effective, yet simple, case.

Given the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(

*t*

_{0}) =

*y*

_{0}

if we want to estimate y(*t*_{1}), we set *h* = *t*_{1} − *t*_{0}. Remember that f( *t*, *y* ) gives the slope at the point (*t*, *y*). Thus, we can find the slope at (*t*_{0}, *y*_{0}):

*K*

_{0}= f(

*t*

_{0},

*y*

_{0})

Next, we use this slope to estimate y(*t*_{0} + *h*/2) ≈ *y*_{0} + ½ *h**K*_{0} and sample the slope at this intermediate point:

*K*

_{1}= f(

*t*

_{0}+ ½

*h*,

*y*

_{0}+ ½

*h*

*K*

_{0})

Using this new slope, we estimate y(*t*_{0} + *h*/2) ≈ *y*_{0} + ½ *h**K*_{1} and sample the slope at this new point:

*K*

_{2}= f(

*t*

_{0}+ ½

*h*,

*y*

_{0}+ ½

*h*

*K*

_{1})

Finally, we use this last approximation of the slope to estimate y(*t*_{1}) = y(*t*_{0} + *h*) ≈ *y*_{0} + *h**K*_{2} and sample the slope at this point:

*K*

_{3}= f(

*t*

_{0}+

*h*,

*y*

_{0}+

*h*

*K*

_{2})

All four of these slopes, *K*_{0}, *K*_{1}, *K*_{2}, and *K*_{3}
approximate the slope of the solution on the interval [*t*_{0}, *t*_{1}],
and therefore we take the following weighted average:

Therefore, we approximate y(*t*_{1}) by

Note how the two values K_{1} and K_{2} both approximate
the slope at the midpoint. Thus, the averaging is very similar to
Simpson's rule which also places a weight of 4 on the midpoint.

# Graphical Representation

To view this graphically, let's consider the IVP

^{(1)}(

*t*) = ½y(

*t*) -

*t*+ 1

y(0) = 0.5

Figure 1 shows Euler's approximation of y(1), that is, following the
slope *K*_{0} = f(*t*_{0}, *y*_{0}) = 1.25, we
get that y(0.5) ≈ 0.5 + 0.5 ⋅ 1.25 = 1.125.

Figure 1. Euler's approximation to the value of y(0.5).

If we sample the slope at (0.5, 1.125), we find that it is
*K*_{1} = f(0.5, 1.125) = 1.0625.
This slope is shown as a light-blue line in Figure 2.

Figure 2. The slope at the point (0.5, 1.125).

We can use this second slope to approximate
y(0.5) ≈ 0.5 + 0.5 ⋅ *K*_{1} = 1.03125. This
is shown in Figure 3.

Figure 3. Following the slope *K*_{1} out a distance *h*/2 = 0.5.

If we sample the slope at (0.5, 1.03125), we find that it is
*K*_{2} = f(0.5, 1.03125) = 1.015625.
This slope is shown as a gold line in Figure 4.

Figure 4. The slope at the point (0.5, 1.03125).

Finally, we can use this third slope to approximate
y(1.0) ≈ 0.5 + 1 ⋅ *K*_{2} = 1.515625. This
is shown in Figure 5.

Figure 5. Following the slope *K*_{2} out a distance *h* = 1.

If we sample the slope at (1, 1.515625), we find that it is
*K*_{3} = f(l, 1.515625) = 0.7578125.
This slope is shown as a magenta line in Figure 4.

Figure 6. The slope at the point (1, 1.515625).

The weighted average of these four slopes is 1.02734375.
Therefore, we approximate y(t_{1}) by following this average slope out
a distance 1: y_{1} = y_{0} + 1.02734375 ⋅ 1 = 1.527343750. This is
shown in Figure 7. Note the significant improvement in the approximation as
compared to that found by Heun's method.

Figure 7. 4th-order Runge Kutta approximation to y(1).

# HOWTO

# Problem

Given the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(t

_{0}) =

*y*

_{0}

approximate y(*t*_{1}).

# Assumptions

The function f(t, y) should be continuous in both variables.

# Tools

We will use Taylor series.

# Initial Conditions

Set *h* = *t*_{1} − *t*_{0}.
Let *y*_{1} be the approximation of y(*t*_{1}).

# Process

Define

*K*_{0} = f( *t*_{0}, *y*_{0} )

*K*_{1} = f( *t*_{0} + ½ *h*, *y*_{0} + ½ *h* *K*_{0} )

*K*_{2} = f( *t*_{0} + ½ *h*, *y*_{0} + ½ *h* *K*_{1} )

*K*_{3} = f( *t*_{0} + *h*, *y*_{0} + *h* *K*_{2} )

Then let

# Examples

# Example 1

Given the IVP y^{(1)}(*t*) = 1 - *t* y(*t*) with
y(0) = 1, approximate y(1) with one step.

First, let *t*_{0} = 0, *y*_{0} = 1, and *h* = 1.
Thus, we calculate

*K*_{0} = f(0, 1) = 1

*K*_{1} = f(0.5, 1.5) = 0.25

*K*_{2} = f(0.5, 1.125) = 0.4375

*K*_{3} = f(1, 1.4375) = -0.4375

(*K*_{0} + 2*K*_{1} + 2*K*_{2} + *K*_{3})/6 = 0.3229166667

Therefore, approximation is *y*_{1} = *y*_{0} + *h* 0.3229166667 = 1.3229166667. The actual value is 1.331309118 and
therefore the absolute error is approximately 0.0084, significantly smaller than the error
using Heun's method.

# Example 2

Given the same IVP shown in Example 1, approximate y(0.5).

*K*_{0} = f(0, 1) = 1

*K*_{1} = f(0.25, 1.25) = 0.6875

*K*_{2} = f(0.25, 1.171875) = 0.70703125

*K*_{3} = f(0.5, 1.353515625) = 0.3232421875

(*K*_{0} + 2*K*_{1} + 2*K*_{2} + *K*_{3})/6 = 0.6853841146

Therefore, approximation is *y*_{1} = *y*_{0} + *h* = 1 + 0.5⋅0.6853851146 = 1.342692057. The actual value is 1.342841185 and
therefore the absolute error is approximately 0.00015, or approximately 1/56th the error
when we used *h* = 1. This is significantly better than the estimated 1/32.

# Example 3

Repeat Examples 1 and 2 but with with the initial value y(0.5) = 2.5 and approximating y(1.5) and y(1.0).

First, let *t*_{0} = 0.5, *y*_{0} = 2.5, and *h* = 1. Therefore

*K*_{0} = f(0.5, 2.5) = -0.25

*K*_{1} = f(1, 2.375) = -1.375

*K*_{2} = f(1, 1.8125) = -0.8125

*K*_{3} = f(0.5, 1.6875) = -0.153125

(*K*_{0} + 2*K*_{1} + 2*K*_{2} + *K*_{3})/6 = -1.026041667

Therefore, our approximation is *y*_{1} = *y*_{0} + *h*(-1.026041667) = 1.473958333. The actual value is 1.502483616 and
therefore the absolute error is approximately 0.0285.

Next, let *t*_{0} = 0.5, *y*_{0} = 2.5, and *h* = 0.5. Therefore

*K*_{0} = f(0.5, 2.5) = -0.25

*K*_{1} = f(1, 2.4375) = -0.828125

*K*_{2} = f(1, 2.29296875) = -0.719726562

*K*_{3} = f(0.5, 2.140136719) = -1.140136719

(*K*_{0} + 2*K*_{1} + 2*K*_{2} + *K*_{3})/6 = -0.7476399739

Therefore, our approximation is *y*_{1} = *y*_{0} + *h*(-0.7476399739) = 2.126180013. The actual value is 2.126611964 and
therefore the absolute error is approximately 0.00043. Again, the error
is approximately 1/66th that of the previous calculation where *h* = 1.

If we tabulate the errors for various values of *h*, as is shown
in Table 1, we note that as *h* gets smaller, the error begins
to drops by a factor 1/32 each time we divide *h* by two.

Table 1. Errors when approximating y(*t*_{0} + *h*) for decreasing values of *h*.

h | Approximation of y(0.5 + h) | Error |
---|---|---|

1 | 1.473958333 | 0.0285 |

0.5 | 2.126180013 | 0.00043 |

0.25 | 2.368600686 | 0.0000089 |

0.125 | 2.450805650 | 0.00000023 |

0.0626 | 2.479808899 | 0.0000000064 |

0.03125 | 2.491036662 | 0.00000000019 |

0.015625 | 2.49580492 | 0.0000000000058 |

You will note that with the last step, the error goes down by approximately a factor of 1/32.7.

# Engineering

To be completed.

# Error

Finding the error of the 4th-order Runge-Kutta method
is beyond the scope of this course. It suffices to state
that it is **O**(*h*^{5}).

# Questions

# Question 1

Given the IVP

^{(1)}(

*t*) = 1 - 0.25 y(

*t*) + 0.2 t

y(0) = 1

approximate y(0.5), y(1), and y(1.5) using 4th-order Runge-Kutta.

Answer: 1.376499430, 1.755761719, 2.137469483.

# Question 2

Given the same ODE as in Question 1, but with the initial condition y(1) = 2, approximate y(1.5) and y(2.0).

Answer: 2.352998861 and 2.711523438

# Matlab

To be completed.

# Maple

Given the IVP y^{(1)}(*t* = *t* y(*t*) - *t*^{2} + 1 with y(0) = 3, to approximate y(1), we can do:

t[0] := 0; t[1] := 1; y[0] := 3; f := (t, y) -> t*y - t^2 + 1; h := t[1] - t[0]; K[0] := f(t[0], y[0] ); K[1] := f(t[0] + h/2, y[0] + h/2*K[0]); K[2] := f(t[0] + h/2, y[0] + h/2*K[1]); K[3] := f(t[0] + h , y[0] + h*K[2]); y[1] := y[0] + h*(K[0] + 2*K[1] + 2*K[2] + K[3])/6;

Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.