# Taha Azzaoui - 2018.03.17

### Introduction

When last we discussed the topic of numerical methods, we touched on Euler’s method for generating a numerical solution to the initial value problem $\frac{dy}{dx} = f(x,y)$ given y(x0) = y0. We concluded by briefly highlighting the tradeoff between simplicity and accuracy that is exhibited in Euler’s method.

One of the problems with Euler’s method is that it fails to effectively consider the curvature of the solution. Since the method essentially uses a line to approximate a curve, it performs quite well when the solution curve is linear. In fact, the estimate turns out to be exact when the solution curve is exactly linear. On the other hand, this also renders Euler’s method ineffective in the context of highly nonlinear solution curves, such as those that arise in oscillating and non-conservative systems.

The following is an outline of an improvement upon Euler’s Method, developed in concert by German Mathematicians Carl Runge and Martin Kutta. Note that although there exist a variety of Runge-Kutta methods, we limit our discussion to the classic fourth order method, which, for the sake of brevity, we refer to as RK4.

### The Problem

Our problem remains the same as before. We seek a solution to the first order ODE: $\frac{dy}{dx} = f(x,y)$ with some known initial value y(x0) = y0. As before, we assume the necessary conditions for the existence of a unique solution, namely that f and fy are continuous.

### The Method

The idea is to take the same iterative approach as we did before, however, this time, we split our interval into four steps, calculating an intermediate value ki at the ith step. After calculating each ki, we calculate our next step as a weighted average of them all. The following is an informal description of each of the intermediate values ki.

• k1: The slope at the current time-step
• k2: The slope at the midpoint of the current time step (halfway between the current one and the next)
• k3: An additional estimate of the midpoint slope, taking into consideration the value of $k_2$
• k4: Our final estimate of the slope at the end of the current time step.

More formally, let ϵ be our step size…

• k1 = f(xi, yi)
• $k_2 = f(x_i + \frac{\epsilon}{2}, y_i + \epsilon \frac{k_1}{2})$
• $k_3 = f(x_i + \frac{\epsilon}{2}, y_i + \epsilon \frac{k_2}{2})$
• k4 = f(xi + ϵ, yi + ϵk3)

Our next step is then a weighted average of the preceding terms …
$$y_{i + 1} = y_i + \frac{\epsilon}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$
At first glance, the formulation of these steps, along with their weighted average seems to be a bit arbitrary. While the derivation of the method is out of the scope of this post, these numbers have been chosen in retrospect of the error analysis.The weights themselves can be derived by solving the set of equations that result from forcing the error terms in the Taylor Series approximation of the solution to cancel out. For a more rigorous discussion on the coefficients and their effect on convergence, see here.

### An Example

Consider the following ODE: $\frac{dy}{dx} = e^{-x}$, y(0) = 1. Although the algorithm is straightforward, approximating solutions by hand can get tedious fairly quickly. We therefore allow the computer to do what it does best and run our example using the following python implementation of RK4. Note that the code used to generate the following examples/figures is available here.

'''
f = f(x,y)
x_0 = initial x
y_0 = initial y
n = number of subdivisions
h = epsilon
'''
def RK4(f,x_0,y_0,n,h):
x = []
y = []
x_i = x_0
y_i = y_0
for i in range(n):
x.append(x_i)
y.append(y_i)
k1 = f(x_i,y_i)
k2 = f(x_i + 0.5 * h, y_i + 0.5 * h * k1)
k3 = f(x_i + 0.5 * h, y_i + 0.5 * h * k2)
k4 = f(x_i + h, y_i + h * k3)
y_i += h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
x_i += h
return x, y

Below is a plot of the output with the parameters n = 50 and ϵ = 0.2. The points in blue are generated by Euler’s method, those in green are generated by RK4, and the purple line is the actual solution, y = 2 − e − x.

### Remarks

We mentioned earlier that there exist a variety of Runge-Kutta methods. In fact, Euler’s method is the simplest member of the class of Runge-Kutta methods, also known as the First Order Runge-Kutta Method. Recall that if the exact solution curve is linear, Euler’s method gives an exact estimate. In the case of RK4, the estimate will be exact if the solution curve is a polynomial of order 4 or less. In general, it can be shown to be the case that if the exact solution is a polynomial of order n or less, it will be estimated exactly by RKn, an nth order Runge-Kutta Method. That being said, as n grows larger it becomes much more computationally expensive to estimate a solution, due to the fact that increasing the order increases the number of Taylor Series terms used, and therefore the number of function evaluations.