Labels

Numerical solutions to differential equations

Differential equations are equations such as $y' = y$, meaning that the derivative of a function is equal to the function itself. This is a very simple example, with an easy solution, but other differential equation do not have a solution that is this simple. Most differential equations have pretty complex solutions, or no analytical solution at all. In the case of $y' = y$, we can see that a solution is $e^x$, because the derivative of $e^x$ is itself.

However, in the case when there is no analytical solution, we can still find a numerical solution. The Runga-Kutta method is a method that provides an numerical solution to a differential equation. The most simple Runga-Kutta method is Euler's method, which we will explain in this post.

Direction fields

Euler's method is based on evaluating the function to find the slope at a point, and then "walking" into the direction of the slope. This can be visualized with the idea of a direction field. If we for example, take $100$ points on a grid, and evaluate the function, we will find the slope at every point. If we then plot the slope, we will get a graph that is called a direction field.

If for example, we have a differential equation $y' = x^2 + y^2 - 1$, then we evaluate the function at a set of points, such as $(0, 0), (0, 1), (1, 0), (1, 1)$, and so on. If we do this for many points on the interval $x \in [-2, 2], y \in [-2, 2]$, we will get the following plot:


If we look at this plot, and take for example point $(-1.5, -1.5)$, we can simply follow the arrows to find the numerical solution to our differential equation $y' = x^2 + y^2 - 1$. This is exactly what Euler's method does. To get a bit more concrete, let's define Euler's method.

Euler's method

If we have a function $F(t)$, and our initial conditions $(x_0, y_0)$ (our starting point), with a step size $h$, then we can approximate a solution to the differential equation with:

$$y_t = y_{t-1} + h\cdot F(x_{t-1}, y_{t-1})$$

Note that we are approximating the function! Also, based on our initial conditions (starting position) we will get a different solution. A differential equation usually has a family of solutions. If we apply this method to our differential equation $F(x, y) = x^2 + y^2 - 1$, with our initial conditions $x_0 = -2$, and $y_0 = -2$, with a step size $h = 0.2$. The first step becomes:

$$ \begin{align} y_1 &= y_0 + h\cdot F(x_0, y_0) \\ y_1 &= -2 + 0.2\cdot\left( (-2)^2 + (-2)^2 - 1\right) \\ y_1 &= 7\end{align}$$

Then for the next step we take $x_1 = x_0 + h$, and then we calculate $y_2$. If we do this $20$ times, we will end up with the following table.

x y
0 -2.0 -2.00000
1 -1.8 -0.60000
2 -1.6 -0.08000
3 -1.4 0.23328
4 -1.2 0.43616
5 -1.0 0.56221
6 -0.8 0.62543
7 -0.6 0.63166
8 -0.4 0.58346
9 -0.2 0.48354
10 -0.0 0.33831
11 0.2 0.16120
12 0.4 -0.02561
13 0.6 -0.19347
14 0.8 -0.31399
15 1.0 -0.36627
16 1.2 -0.33944
17 1.4 -0.22840
18 1.6 -0.02596
19 1.8 0.28617
20 2.0 0.75055

Now by looking at these values we are not getting any wiser, so here is a graph of the direction field, with these values plotted in.


So there you have it, a method to find a numerical solution to a differential equation. Note that Euler's method is the simplest Runge-Kutta method. If greater accuracy is required, it is advised to use a fourth order Runge-Kutta method (RK4), also knowns as the classic Runge-Kutta method.

Code used for plotting

For completeness sake, and for myself to find this back, this is the Python code that is used to create the graphs. It requires the following modules: numpy, matplotlib, pandas.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

fig = plt.figure(figsize=(7,7))

# Vector field
X, Y = np.meshgrid(np.linspace(-2, 2, 20), np.linspace(-2, 2, 20))
U = 1.0
V = X*X + Y*Y - 1
# Normalize arrows
N = np.sqrt(U ** 2 + V ** 2)
U = U / N
V = V / N
plt.quiver(X, Y, U, V, angles="xy")

vf = lambda x, t: (1, x[0] + x[1] - 1)

t = np.linspace(t0, tEnd, 100)
for y0 in np.linspace(-5.0, 0.0, 10):
    y_initial = [y0, -10.0]
    y = odeint(vf, y_initial, t)
    plt.plot(y[:, 0], y[:, 1], "-")

#plt.grid(alpha=0.7, ls='dashed')
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")

# Euler's method for find a numerical solution:
# Based on the formula:
#    y_n = y_n+1 + h * F(x_n-1, y_n-1), and x_n = x_n-1 + h
#    with starting conditions x_0, y_0 and step size h.

# Initial conditions
xn = -2
yn = -2
h = 0.2 # Step size
xs = [yn]
ys = [xn]

# Function f(x, y)
f = lambda x, y: x*x + y*y - 1

# Perform Euler's method
n = 20
for _ in range(n):
    yn = yn + h*f(xn, yn)
    xn += h
    xs.append(xn)
    ys.append(yn)
    
# Plot the solution
plt.title("Direction field of $x^2 + y^2 - 1$")
plt.plot(xs, ys, lw=3, c='red')

# Show table
df = pd.DataFrame(data={'x': xs, 'y': ys})
df.round(5)
Labels: