Ordinary Differential Equation Solvers: Runge-Kutta Methods

Christina Lee

Category: Numerics

So what's an Ordinary Differential Equation?

Differential Equation means we have some equation (or equations) that have derivatives in them.

The ordinary part differentiates them from partial differential equations (the ones with curly $\partial$ derivatives). Here, we only have one independent variable (let's call it $t$), and one or more dependent variables (let's call them $x_1, x_2, ...$). In partial differential equations, we can have more than one independent variable.

This ODE can either be written as a system of the form $$ \frac{d x_1}{dt}=f_1(t,x_1,x_2,...,x_n) $$ $$ \frac{d x_2}{dt}=f_2(t,x_1,x_2,...,x_n) $$ ... $$ \frac{d x_n}{dt}=f_n(t,x_1,x_n,...,x_n) $$ or a single n'th order ODE of the form $$ f_n(t,x) \frac{d^n x}{dt^n}+...+f_1(t,x) \frac{dx}{dt}+f_0(t,x)=0 $$ that can be rewritten in terms of a system of first order equations by performing variable substitutions such as $$ \frac{d^i x}{dt^i}=x_i $$

Though STEM students such as I have probably spent thousands of hours pouring of ways to analytically solve both ordinary and partial differential equations, unfortunately, the real world is rarely so kind as to provide us with an exactly solvable differential equation. At least for interesting problems.

We can sometimes approximate the real world as an exactly solvable situation, but for the situation we are interested in, we have to turn to numerics. I'm not saying those thousand different analytic methods are for nothing. We need an idea ahead of time of what the differential equation should be doing, to tell if it's misbehaving or not. We can't just blindly plug and chug.

Today will be about introducing four different methods based on Taylor expansion to a specific order, also known as Runge-Kutta Methods. We can improve these methods with adaptive stepsize control, but that is a topic for another time, just like the other modern types of solvers such as Richardson extrapolation and predictor-corrector.

Nonetheless, to work with ANY computational differential equation solver, you need to understand the fundamentals of routines like Euler and Runge-Kutta, their error propagation, and where they can go wrong. Otherwise, you might misinterpret the results of more advanced methods.

WARNING: If you need to solve a troublesome differential equation for a research problem, use a package, like DifferentialEquations. These packages have much better error handling and optimization.

Let's first add our plotting package and colors.

We will be benchmarking our solvers on one of the simplest and most common ODE's,

\begin{equation} \frac{d}{d t}x=x \;\;\;\;\;\;\; x(t)=C e^t \end{equation}

Though this only has one dependent variable, we want to structure our code so that we can accommodate a series of dependent variables, $y_1,y_2,...,y_n$, and their associated derivative functions. Therefore, we create a function for each dependent variable, and then push it into an array declared as type Function.

Euler's Method

First published in Euler's Instutionum calculi integralis in 1768, this method gets a lot of milage, and if you are to understand anything, this method is it.

We march along with step size $h$, and at each new point, calculate the slope. The slope gives us our new direction to travel for the next $h$.

We can determine the error from the Taylor expansion of the function $$ x_{n+1}=x_n+h f(x_n,t) + \mathcal{O}(h^2). $$ In case you haven't seen it before, the notation $\mathcal{O}(x)$ stands for "errors of the order x". Summing over the entire interval, we accumuluate error according to $$N\mathcal{O}(h^2)= \frac{x_f-x_0}{h}\mathcal{O}(h^2)=h, $$ making this a first order method. Generally, if a technique is $n$th order in the Taylor expansion for one step, its $(n-1)$th order over the interval.

Implicit Method or Backward Euler

If $f(t,x)$ has a form that is invertible, we can form a specific expression for the next step. For example, if we use our exponential, \begin{equation} x_{n+1}=x_n+ h f(t_{n+1},x_{n+1}) \end{equation} \begin{equation} x_{n+1}-h x_{n+1}=x_n \end{equation} \begin{equation} x_{n+1}=\frac{x_n}{1-h} \end{equation}

This expression varies for each differential equation and only exists if the function is invertible.

2nd Order Runge-Kutta

So in the Euler Method, we could just make more, tinier steps to achieve more precise results. Here, we make bettter steps. Each step itself takes more work than a step in the first order methods, but we win by having to perform fewer steps.

This time, we are going to work with the Taylor expansion up to second order, \begin{equation} x_{n+1}=x_n+h f(t_n,x_n) + \frac{h^2}{2} f^{\prime}(t_n,x_n)+ \mathcal{O} (h^3). \end{equation}

Define \begin{equation} k_1=f(t_n,x_n), \end{equation} so that we can write down the derivative of our $f$ function, and the second derivative (curvature), of our solution, \begin{equation} f^{\prime}(t_n,x_n)=\frac{f(t_n+h/2,x_n+h k_1/2)-k_1}{h/2}+\mathcal{O}(h^2). \end{equation} Plugging this expression back into our Taylor expanion, we get a new expression for $x_{n+1}$ \begin{equation} x_{n+1}=x_n+hf(t_n+h/2,x_n+h k_1/2)+\mathcal{O}(h^3) \end{equation}

We can also interpret this technique as using the slope at the center of the interval, instead of the start.

4th Order Runge-Kutta

Wait! Where's 3rd order? There exists a 3rd order method, but I only just heard about it while fact-checking for this post. RK4 is your dependable, multi-purpose workhorse, so we are going to skip right to it.

$$ k_1= f(t_n,x_n) $$$$ k_2= f(t_n+h/2,x_n+h k_1/2) $$$$ k_3 = f(t_n+h/2, x_n+h k_2/2) $$$$ k_4 = f(t_n+h,x_n+h k_3) $$$$ x_{n+1}=x_n+\frac{h}{6}\left(k_1+2 k_2+ 2k_3 +k_4 \right) $$

I'm not going to prove here that the method is fourth order, but we will see numerically that it is.

Note: I premultiply the $h$ in my code to reduce the number of times I have to multiply $h$.

This next function merely iterates over a certain number of steps for a given method.

Scaling of the Error

I talked above about the error scaling either as $h,h^2$, or $h^4$. I won't just talk but here will numerically demonstrate the relationship as well. For a variety of different step sizes, the below code calculates the final error for each method. Then we will plot the error and see how it scales.

Arbitrary Order

While I have presented 4 concrete examples, many more exist. For any choice of variables $a_i, \beta_{i,j},a_i$ that fulfill $$ x_{n+1}=x_n+h\left(\sum_{i=1}^s a_i k_i \right)+ \mathcal{O}(h^p) $$ with $$ k_i=f\left(t_n+\alpha_i h,x_n+h\left(\sum_{j=1}^s \beta_{ij} k_j \right) \right) $$ we have a Runge-Kutta method of order $p$, where $p\geq s$. The Butcher tableau provides a set of consistent coefficients.

Stay tuned for when we tuned these routines to the stiff van der Pol equations!