Skip to main content

Generalized Explicit Runge-Kutta

Project description

Generalized Explicit Runge-Kutta (GERK)

A package for the curious mathematicians and engineers who want to experiment the Runge-Kutta method with their own coefficients.

PyPI link

GitHub link

Motivation

I decided to create this package to give mathematicians or engineers the freedom to create and use the Runge-Kutta method of their choice, as opposed to be being locked in to whatever a package provides you.

Runge-Kutta Overview

Runge-Kutta methods aim to numerically solve ordinary differential equations of the form:

$$ \frac{\text{d}y}{\text{d}x}=f(x,y) $$

with a given initial condition $(x_0, y_0)$.

We define a member of the Runge-Kutta family with a Butcher tableau:

genrk

The above tableau is often abbreviated to

abbrk

The Butcher tableau presents the categories of coefficients that will be used in our Runge-Kutta method.

The $n^\text{th}$ evaluation of the solution will be denoted as $(x_n, y_n)$. We also define $h$ as the time step i.e. the step size from the previous approximation to the next, and therefore

$$ x_{n+1} = x_{n} + h $$

Before defining $y_n$, we must familiarize ourselves with the array $k$. We define the $i^{\text{th}}$ row of $k$ at $(x_n, y_n)$ as:

$$ k_i(x_n, y_n) = f\left(x_n + c_i h, y_n + \sum_{j=1}^{i-1} a_{ij}k_j(x_n, y_n)\right) $$

where $f$ is the function defined in the differential equation above. Note the recursion in the second argument of $f$ where we sum rows preceding $k_{i}$ and apply a scale factor of $a_{ij}$.

We now have everything we need to calculate $y_{n+1}$:

$$ y_{n+1} = y_{n} + h \sum_{i=1}^{r}b_i k_{i}(x_n, y_n) $$

The idea is to calculate various slopes at point $y_n$ to ascertain a weighted of the ascent (or descent) and add it to the previous approximation.

What is Gerk?

Most packages for the Runge-Kutta method usually have the coefficients $a_{ij}$, $b_i$ and $c_i$ determined beforehand for known methods such as the Forward-Euler method, the 1/4 rule, the 3/8 rule etc, but do not allow one to customerize their own Runge-Kutta.

Gerk is an easy interface to allow the user to determine their own coefficient values for the Runge-Kutta method. The package can return the approximations in the form of an array, a plot of the approximated curve (using matplotlib) and also compare with real values (if available) to produce errors and even an efficiency graph!

How to use Gerk

We can simply import Gerk in the following way:

from gerk import Gerk

Parameters

As seen in mathematics above, there is quite a bit information required to execute the Runge-Kutta method. This has been broken down into arguments to be passed into the Gerk class:

  • a The $A$ matrix in the Butcher tableau. This must be a lower triangular matrix that is formatted as a list of lists that contain floats, or integers
  • b The $b$ array. Must be a list of floats, decimals or integers
  • c The $c$ array. Must be a list of floats, decimals or integers
  • initial A tuple that acts as the coordinate of the initial condition values $(x_0, y_0)$
  • terminal The value of $x$ for which we terminate the Runge-Kutta method
  • timesteps The number of times steps you want to apply on the from the starting point $x_0$ to final. Must be an integer
  • func The function expression to be numerically integrated
  • enforce_rules (Optional) A boolean to enforce conventional Runge-Kutta rules

Conditions

There is no consensus to what conditions must hold regarding the coefficients you choose for your method, however, some known Runge-Kutta methods do consistently conform to some known conditions.

These conditions are

$$\sum_{i=1}^{r}b_i=1 \ \ \ \ \sum_{i=1}^{r}b_ic_i = 1/2 \ \ \ \ \sum_{j=1}^{r}a_{ij} = c_i$$

which can be enforced by the enforce_rules respectively. These parameters have been set to False by default to allow the user to freely explore and experiment.

Example

Let us numerically solve the following initial value problem:

$$ \frac{\text{d}y}{\text{d}x} = y $$

with the initial value $(0, 1)$. We want to apply a Runge-Kutta method with the following Butcher tableau:

The $A$ lower triangular matrix in the Butcher tableau above can be implemented in the following way:

a = [ 
        [1/3],
        [-1/3, 1],
        [1, -1, 1]
]

Now for the b and c vectors:

b = [1/8, 3/8, 3/8, 1/8]

c = [0, 1/3, 2/3, 1]

Finally, we can define the function with a lambda function:

func = lambda x, y: y

Now we are ready to run gerk() and plot its output:

import matplotlib.pyplot as plt
from gerk import gerk

x, y = gerk(
        a=a, 
        b=b, 
        c=c, 
        initial=(0, 1), 
        timesteps=10000, 
        terminal=1,
        func=lambda x, y: y,
    )

plt.plot(x, y, color="r")
plt.show()
gerk_1

Adaptive Runge-Kutta Methods

There is an alternate way to utilise the Runge-Kutta method by employing an additional distinct $b$ array. The Butcher tableau for such methods take the form:

adaptive_butcher

where $b^*_i$ is the additional $b$ array.

This method is not too disimilar to the original Runge-Kutta method. We calculate the $k$'s in the same way as before, but there is an extra step when evaluating $y_{n+1}$.

Although we do calculate $y_{n+1}$ in same way outlined above, we also calculate $\hat{y}_{n+1}$:

$$ y_{n+1} = y_{n} + h \sum_{i=1}^{r}b_{1i}\cdot k_{i}(x_n, y_n) \ \ \ \ \ \ \ \ \ \hat{y_{n+1}} = y_{n} + h \sum_{i=1}^{r} b_{2i}\cdot k_{i}(x_n, y_n) $$

Note that the calculation for $\hat{y}$ requires the use of $y_n$ and not $\hat{y}_n$.

For every time step, we calculate the following error:

$$ E := \left|y_{n+1}-\hat{y}_{n+1}\right| $$

Now we define the tolerance, $\mathcal{E}$, which will act as the maximum acceptable value for $E$.

If $E<\mathcal{E}$, then we accept the value of $y_{n+1}$ and we increment $x_{n}$ by $h$ and start the process again for $x_{n+1}$ and $y_{n+1}$ as normal.

However, if $E\geq\mathcal{E}$, we will need to adjust the value of $h$ and redo the calculation with this renewed time step value in an attempt to satisfy the condition $E<\mathcal{E}$.

The value of $h$ will be adjusted as follows:

$$ h \rightarrow 0.9\cdot h \cdot\sqrt[n]{\frac{h}{\mathcal{E}}} $$

where $n=\min\left(p,q\right)$, where $p$ and $q$ are the orders $^1$ of $b_i$ and $b^*_i$ respectively.

$^1$ A Runge-Kutta method with matrix $A$ and arrays $b$ and $c$ has order $p$ if

$$ \sum_{i=1}^{p}b_i = 1 \ \ \ \ \sum_{i=1}^{p}b_ic_i = 1/2 \ \ \ \ \sum_{j=1}^{p}a_{ij} = c_i \ \ \ \ $$

Adaptive Runge-Kutta example

In this example, we will try to numerically integrate

$$ \frac{\text{d}y}{\text{d}x}=-2xy $$

with initial conditions $(-1, e^{-1})$. Note that the exact solution is $y=e^{-x^2}$.

Here we will use the Bogacki–Shampine (BS23) method which has the following Butcher tableau:

adapt_butcher

For the adaptive Runge-Kutta method, we will use adaptive_gerk(). This method's paramters vary slightly from gerk():

  • a The $A$ matrix in the Butcher tableau. This must be a lower triangular matrix that is formatted as a list of lists that contain floats, or integers
  • b_1 The $b_1$ array. Must be a list of floats, decimals or integers
  • b_2 The $b_2$ array. Must be a list of floats, decimals or integers
  • c The $c$ array. Must be a list of floats, decimals or integers
  • initial A tuple that acts as the coordinate of the initial condition values $(x_0, y_0)$
  • terminal The value of $x$ for which we terminate the Runge-Kutta method
  • timesteps The number of times steps you want to apply on the from the starting point $x_0$ to final. Must be an integer
  • func The function expression to be numerically integrated
  • enforce_rules (Optional) A boolean to enforce conventional Runge-Kutta rules. Defaulted to False
  • tolerance (Optional) A float representing the maximum threshold of the adaptive step. Defaulted to 1e-4
from math import exp
import matplotlib.pyplot as plt
from gerk import adaptive_gerk

a = [
    [1/3],
    [-1/3, 1],
    [1, -1, 1]
]
b_1 = [2/9, 1/3, 4/9, 0]
b_2 = [7/24, 1/4, 1/3, 1/8]
c = [0, 1/2, 3/4, 1]

x, y = adaptive_gerk(
        a=a, 
        b_1=b_1, 
        b_2=b_2, 
        c=c, 
        initial=(-1, exp(-1)), 
        timesteps=10000, 
        terminal=1,
        func=lambda x, y: -2 * x * y,
    )

plt.plot(x, y, color="r")
plt.show()
gerk_2

In this case, the Runge-Kutta approximation is so accurate that we can barely see the exact curve!

Efficiency graphs are not available for the adaptive Runge-Kutta method as the number discretizations is not constant.

Future work

  • To further generalise the Runge-Kutta method with an indefinite number of $b$ and $c$ arrays

Project details


Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Source Distribution

gerk-1.0.1.tar.gz (6.9 kB view details)

Uploaded Source

Built Distribution

gerk-1.0.1-py3-none-any.whl (7.7 kB view details)

Uploaded Python 3

File details

Details for the file gerk-1.0.1.tar.gz.

File metadata

  • Download URL: gerk-1.0.1.tar.gz
  • Upload date:
  • Size: 6.9 kB
  • Tags: Source
  • Uploaded using Trusted Publishing? No
  • Uploaded via: twine/5.1.1 CPython/3.10.6

File hashes

Hashes for gerk-1.0.1.tar.gz
Algorithm Hash digest
SHA256 f53f3d7a569b65d176da51ee48ab5bd29dbabe74b296668ca44397fed1f9b315
MD5 c846e24889864e7b33c04bb490ab5c5c
BLAKE2b-256 2b9393be985f673069f56ff2e5611c627afe7736af7bb1181b42262047335ad8

See more details on using hashes here.

File details

Details for the file gerk-1.0.1-py3-none-any.whl.

File metadata

  • Download URL: gerk-1.0.1-py3-none-any.whl
  • Upload date:
  • Size: 7.7 kB
  • Tags: Python 3
  • Uploaded using Trusted Publishing? No
  • Uploaded via: twine/5.1.1 CPython/3.10.6

File hashes

Hashes for gerk-1.0.1-py3-none-any.whl
Algorithm Hash digest
SHA256 0ad38e4ab5f74e37e495b99a05ded90235bfe27b9beed34f66feb8a12fb68109
MD5 8ef7f5afe0881395094cce718a5dcaba
BLAKE2b-256 1eeb6bbfebb215f0635c39bf208a24c662561ea23360bba7f05c28302b394c1b

See more details on using hashes here.

Supported by

AWS AWS Cloud computing and Security Sponsor Datadog Datadog Monitoring Fastly Fastly CDN Google Google Download Analytics Microsoft Microsoft PSF Sponsor Pingdom Pingdom Monitoring Sentry Sentry Error logging StatusPage StatusPage Status page