Up to now, we have seen three techniques for reducing the error for solving an initial value problem of the form:

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(

*a*) =

*y*

_{0}

where we estimated for y(*b*) with a single step (or iteration).
To get a better answer, use the same strategy as we did for the
composite trapezoidal rule by breaking the interval up into smaller
sub intervals and applying either Euler's, Heun's or the 4th-order
Runge Kutta methods on each subinterval.

# Background

Useful background for this topic includes:

- 4. Linear Algebra
- 14.1 Euler's Method
- 14.2 Heun's Method
- 14.3 4th-Order Runge-Kutta Method

# References

See the corresponding references in the last three background topics.

# Theory

Given the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(

*a*) =

*y*

_{0}

where we want to estimated y(*b*) for *b* > *a*,
we can apply the same strategy we used for the composite-trapezoidal rule:

Break the interval [*a*, *b*] into *n* sub-intervals
and define *h* = (*b* - *a*)/*n*. Then set
*t*_{i} = *a* + *ih* for *i* = 0, 1, ..., *n*.
Therefore *t*_{0} = *a* and *t*_{n} = *b*.

Now, *y*_{0} is the initial value so let
*y*_{i} represent the approximation of
y(*t*_{i}) for *i* = 1, 2, ..., *n*.

Using any of the three techniques we've seen, Euler's, Heun's, or 4th-order Runge Kutta, we may now proceed as follows:

For *i* from 1 to *n*:

Approximate *y*_{i} using *y*_{i − 1}.

More specifically:

### Multiple-step Euler's Method

Calculate

*y*

_{i}=

*y*

_{i − 1}+

*h*f(

*t*

_{i − 1},

*y*

_{i − 1})

for *i* = 1, 2, ..., *n*.

### Multiple-step Heun's Method

Calculate

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

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

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + *K*_{1})/2.

### 4th-order Runge Kutta

Calculate

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

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

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

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

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + 2 *K*_{1} + 2 *K*_{2} + *K*_{3})/6.

# HOWTO

# Problem

Given the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(

*a*) =

*y*

_{0}

approximate y(*b*).

# Assumptions

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

# Tools

We will use Taylor series and iteration.

# Initial Conditions

Choose a value of *n* ≥ 1. Set *h* = (*b* − *a*)/*n*.

Set *t*_{i} = *a* + *ih* for *i* = 0, 1, ..., *n*
and let *y*_{i} be the approximation of y(*t*_{i})
for *i* = 1, ..., *n*.

# Process

### Multiple-step Euler's Method

For *i* = 1, 2, ..., *n* set

*y*

_{i}=

*y*

_{i − 1}+

*h*f(

*t*

_{i − 1},

*y*

_{i − 1}).

### Multiple-step Heun's Method

For *i* = 1, 2, ..., *n* calculate

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

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

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + *K*_{1})/2.

### 4th-order Runge Kutta

For *i* = 1, 2, ..., *n* calculate

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

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

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

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

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + 2 *K*_{1} + 2 *K*_{2} + *K*_{3})/6.

### Approximating Intermediate Values

If a value of y(*t*) is required for *a* < *t* < *b*
where *t* ≠ *t*_{i} for
any *i*, choose the surrounding four points (three if *t* < *t*_{1}
or *t* > *t*_{n − 1}), find the interpolating
polynomial and evaluate this polynomial at the point *t*. Alternatively, calculate
the appropriate cubic spline where the derivatives at the end points are given
by f(*a*, *y*_{0}) and f(*b*, *y*_{n}),
respectively.

# Examples

# Example 1

Perform four steps of each of Euler's method, Heun's method, and 4th-order Runge Kutta on the IVP

y^{(1)}(*t*) = y(*t*) t + t - 1

y(0) = 1

In this case, the function f(*t*, *y*) = *y* *t* + *t* - 1.
Let *h* = (1 - 0)/4 = 0.25 and thus *t*_{0} = 0, *t*_{1} = 0.25, *t*_{2} = 0.5,
*t*_{3} = 0.75, *t*_{4} = 1 and *y*_{0} = 1. Then

### Euler's method

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

*y*_{2} = *y*_{1} + *h* f(*t*_{1}, *y*_{1}) = 0.609375

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

*y*_{4} = *y*_{3} + *h* f(*t*_{3}, *y*_{3}) = 0.603149414

### Heun's method

K_{0} = f(*t*_{0}, *y*_{0}) = -1

K_{1} = f(*t*_{0} + *h*, *y*_{0} + *h* K_{0}) = -0.5625000000

y_{1} = *y*_{0} + ½ *h* (K_{0} + K_{1}) = 0.8046875000

K_{0} = f(*t*_{1}, *y*_{1}) = -0.5488281250

K_{1} = f(*t*_{1} + *h*, *y*_{1} + *h* K_{0}) = -0.1662597656

y_{2} = *y*_{1} + ½ *h* (K_{0} + K_{1}) = 0.7153015137

K_{0} = f(*t*_{2}, *y*_{2}) = -0.1423492432

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

y_{3} = *y*_{2} + ½ *h* (K_{0} + K_{1}) = 0.7299810648

K_{0} = f(*t*_{3}, *y*_{3}) = 0.297485799

K_{1} = f(*t*_{3} + *h*, *y*_{3} + *h* K_{3}) = 0.804352515

y_{4} = *y*_{3} + ½ *h* (K_{0} + K_{1}) = 0.8677108540

### 4th-order Runge Kutta

K_{0} = f(*t*_{0}, *y*_{0}) = -1

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

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

K_{3} = f(*t*_{0} + *h*, *y*_{0} + *h*K_{0}) = -0.5476226806

y_{1} = *y*_{0} + *h* (K_{0} + 2 K_{1} + 2 K_{2} + K_{3})/6 = 0.8082167307

K_{0} = f(*t*_{1}, *y*_{1}) = -0.5479458173

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

K_{2} = f(*t*_{1} + ½*h*, *y*_{1} + ½*h*K_{1}) = -0.3382126488

K_{3} = f(*t*_{1} + *h*, *y*_{1} + *h*K_{0}) = -0.1381682158

y_{2} = *y*_{1} + *h* (K_{0} + 2 K_{1} + 2 K_{2} + K_{3})/6 = 0.7224772847

K_{0} = f(*t*_{2}, *y*_{2}) = -0.1387613576

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

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

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

y_{3} = *y*_{2} + *h* (K_{0} + 2 K_{1} + 2 K_{2} + K_{3})/6 = 0.7417768882

K_{0} = f(*t*_{3}, *y*_{3}) = 0.306332666

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

K_{2} = f(*t*_{3} + ½*h*, *y*_{3} + ½*h*K_{1}) = 0.585037893

K_{3} = f(*t*_{3} + *h*, *y*_{3} + *h*K_{0}) = 0.888036361

y_{4} = *y*_{3} + *h* (K_{0} + 2 K_{1} + 2 K_{2} + K_{3})/6 = 0.8867587481

The correct answer is y(1) = 0.8867564079.

The plot of the solution and the field plot are shown in Figure 1.

Figure 1. The field plot and the solution.

The approximations are shown in Figure 2.

Figure 2. The solution and the three approximations: Euler (blue), Heun (light blue), and 4th-order Runge Kutta (red)

# Engineering

To be completed.

# Error

Consider the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*))

y(

*a*) =

*y*

_{0}

where we are trying to approximate y(*b*).
If we apply Euler's method with *n* steps, then the error on the
*i*th interval
is given by ½y^{(2)}(τ_{i})*h*^{2}
where τ_{i} ∈ [*t*_{i − 1}, *t*_{i}]. Thus, the
cumulative error is:

If we factor out a ½*h* and rewrite the other *h* as
(*b* - *a*)/*n*, then we get:

Using the definition of the average value, this simplifies to:

Thus, we can easily generalize this result to see that multiple-step methods
will reduce the order of the error term by a factor
of *h*. Thus,
while a single step of 4th-order Runge Kutta is O(*h*^{5}),
multiple applications of 4th-order Runge Kutta is only O(*h*^{4}).
(Hence the name, *4th-order Runge Kutta*.)
A summary is shown in Table 1.

Table 1. Comparison of rates of convergence for single step and multiple steps.

Method | Single Step | Multiple Steps |
---|---|---|

Euler's method | O(h^{2}) | O(h) |

Heun's method | O(h^{3}) | O(h^{2}) |

4th-order Runge Kutta | O(h^{5}) | O(h^{4}) |

This suggests very strongly that it is that much more important to use the 4th-order Runge Kutta method when using multiple step methods.

To justify this statement, let us look at how the error decreases
as the size of the interval is increased. We cannot, however, simply divide
the interval into *n* sub-intervals and apply each method because
Euler's method requires only one function evaluation, and 4th-order Runge Kutta
requires four. Thus, Tables 1, 2, 3, and 4 look at how the error is reduced with
a constant number of function evaluations.

Table 1. Four function evaluations.

Method | Steps | Error Reduction |
---|---|---|

Euler's | 4 | 1/4 = 0.25 |

Heun's | 2 | 1/2^{2} = 0.25 |

4th-order Runge Kutta | 1 | 1 |

Table 2. Eight function evaluations.

Method | Steps | Error Reduction |
---|---|---|

Euler's | 8 | 1/8 = 0.125 |

Heun's | 4 | 1/4^{2} = 0.0625 |

4th-order Runge Kutta | 2 | 1/2^{4} = 0.0625 |

Table 3. Sixteen function evaluations.

Method | Steps | Error Reduction |
---|---|---|

Euler's | 16 | 1/16 = 0.0625 |

Heun's | 8 | 1/8^{2} = 0.0156 |

4th-order Runge Kutta | 4 | 1/4^{4} = 0.00391 |

Table 4. Thirty-two function evaluations.

Method | Steps | Error Reduction |
---|---|---|

Euler's | 32 | 1/32 = 0.03125 |

Heun's | 16 | 1/16^{2} = 0.00391 |

4th-order Runge Kutta | 8 | 1/8^{4} = 0.000244 |

Thus, while Euler's method with 32 sub-intervals and 4th-order Runge Kutta with 8 intervals both have the same number of evaluations, the error with Euler's method is reduced by only 0.03125 while the error with 4th-order Runge Kutta is reduced by 0.000244.

To visualize this, consider:

y^{(1)}(*t*) = 1 + *y*^{2}(*t*)

y(0) = 0

The correct answer is y(1) = tan(1) = 1.557407725.

With Euler's method, using four steps, we have:

*y*_{1} = 0.25

*y*_{2} = 0.515625

*y*_{3} = 0.8320922852

*y*_{4} = 1.255186678

The absolute error of this answer is 0.30.

With one step of 4th-order Runge Kutta, we have:

*y*_{1} = 1.535847982

The absolute error here is 0.022.

If however, we use 32 function evaluations, we have *y*_{32} = 1.497473902
which has an absolute error of 0.060 which is approximately 1/5 the
error with 4 intervals.

With 8 steps of 4th-order Runge Kutta, we get that *y*_{8} = 1.557402847
which has an absolute error of 0.0000049 which is approximately 0.00022 that
of the approximation with one step.

# Questions

# Question 1

Given the IVP

^{(1)}(

*t*) = 1 - 0.25 y(

*t*) + 0.2 t

y(0) = 1

approximate y(1) using four steps of Euler's method, Heun's method, and 4th-order Runge Kutta.

Answer: 1.754495239, 1.755786729, 1.755760161

# Question 2

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

Answer: 2.708990479, 2.711573457, 2.711520324

# Matlab

To be completed later.

# Maple

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

### Euler's Method

h := (b - a)/n; for i from 0 to n do t[i] := a + i*h; end; y[0] := 3; f := (t, y) -> t*y - t^2 + 1; for i from 1 to n do y[i] := y[i - 1] + h*f(t[i - 1], y[i - 1]); end do; y[n]; # approximation

### Heun's Method

h := (b - a)/n; for i from 0 to n do t[i] := a + i*h; end; y[0] := 3; f := (t, y) -> t*y - t^2 + 1; for i from 1 to n do K[0] := f(t[i - 1], y[i - 1]); K[1] := f(t[i - 1] + h, y[i - 1] + h*K[0]); y[i] := y[i - 1] + h*(K[0] + K[1])/2; end do; y[n]; # approximation

### 4th-order Runge Kutta

h := (b - a)/n; for i from 0 to n do t[i] := a + i*h; end; y[0] := 3; f := (t, y) -> t*y - t^2 + 1; for i from 1 to n do K[0] := f(t[i - 1], y[i - 1]); K[1] := f(t[i - 1] + h/2, y[i - 1] + h/2*K[0]); K[2] := f(t[i - 1] + h/2, y[i - 1] + h/2*K[1]); K[3] := f(t[i - 1] + h, y[i - 1] + h*K[2]); y[i] := y[i - 1] + h*(K[0] + 2*K[1] + 2*K[2] + K[3])/6; end do; y[n]; # approximation

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