PyLie

A Python package for solving ordinary differential equations evolving on non-linear manifolds.

This package is distributed with the Python package index. To install it, use

$ pip install pylie

In order to solve an ODE, the differential must first be described in its canonical Lie form – that is, as a mapping from the manifold to the corresponding Lie algebra. For examples, please see below.

Example: Equation evolving on the unit sphere

The complete code is listed at the bottom of this section if you want to copy-paste it, including a definition of A(t, y).

The unit sphere has Lie algebra so(3), consisting of 3-by-3 skew-symmetric matrices (i.e. matrices which satisfy the equation transpose(A) = -A). Ordinary differential equations where the solution space is the unit sphere may be formulated in the form

dy / dt = A(t, y) · y

where A is a skew-symmetric matrix. In order to solve the above equation, you must define the function

def A(t, y):
    # return a 3-by-3 skew-symmetric matrix of type np.ndarray

For instance:

import numpy as np


def A(t, y):
    return np.array(
            [
                [0,                  t,           -0.4 * np.cos(t)],
                [-t,                 0,                0.1 * t    ],
                [0.4 * np.cos(t), -0.1 * t,                0      ]
            ]
        )

You must also decide which numerical scheme you would like to use to solve the equation. Higher-order methods provide a more accurate solution, but are more computationally expensive. For a list of available methods, see available numerical schemes. In this example, we will use the Lie group method corresponding to the fourth order Runge-Kutta method.

To solve the problem, we use the following code:

import numpy as np
import pylie

### Code defining or importing A(t, y) ###

y0 = [0.0, 0.0, 1.0]
t_start = 0
t_end = 5
step_length = 0.1
manifold = "hmnsphere"
method = "RKMK4"
solution = pylie.solve(A, y0, t_start, t_end, step_length, manifold, method)

The variable solution is now a Flow object with two attributes: T, a one-dimensional numpy array containing the times at which the solution is computed, Y, a 3-by-n numpy array where column Y[i, :] is the solution at time T[i]. It is also possible to use indexing directly on the object: solution[i, j] is equivalent to solution.Y[i, j]. If you wish you may also extract the variables Y and T directly by using

# If solution is not yet computed:
Y, T = pylie.solve(A, y0, t_start, t_end, step_length, manifold, method)

# Or, if you followed the example above
Y, T = solution

The following is a suggestion in order to plot the solution:

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.plot(solution[0, :], solution[1, :], solution[2, :])
plt.show()

Full example

This file is also avaiable in /docs/unit_sphere_example.py.

from numpy.testing import assert_almost_equal
import pylie
import numpy as np
import matplotlib.pyplot as plt


def A(t, y):
    return np.array(
        [[0, t, -0.4 * np.cos(t)], [-t, 0, 0.1 * t], [0.4 * np.cos(t), -0.1 * t, 0]]
    )


if __name__ == "__main__":
    y0 = [0.0, 0.0, 1.0]
    t_start = 0
    t_end = 5
    step_length = 0.01
    manifold = "hmnsphere"
    method = "RKMK4"
    solution = pylie.solve(A, y0, t_start, t_end, step_length, manifold, method)

    # Verify that the solution is indeed on the unit sphere
    solution_norm = [np.linalg.norm(solution[:, i]) for i in range(len(solution.T))]
    for val in solution_norm:
        assert_almost_equal(val, 1.0)
    print("Passed test, plotting ...")

    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")
    ax.plot(solution[0, :], solution[1, :], solution[2, :])
    plt.show()

Example: The heavy top

A more detailed explanation is in progress. Until then, the code for the full example is avaible below and in the file /docs/heavy_top_example.py.

from numpy.testing import assert_almost_equal
import pylie
import numpy as np
import matplotlib.pyplot as plt


def heavy_top(
    t, y, principal_moments=np.array([2, 2, 1]), m=1, g=1, chi=np.array([0, 0, 1])
):
    """A formulation of the problem exploiting the Lie-group structure"""
    mu, beta = np.split(y, 2)
    mu_dot = -mu / principal_moments
    beta_dot = -m * g * chi
    return np.hstack((mu_dot, beta_dot))


if __name__ == "__main__":
    y0 = np.array([np.sin(1.1), 0, np.cos(1.1), 1, 0.2, 3])
    t_start = 0
    t_end = 5
    step_length = 0.01
    manifold = "heavytop"
    method = "RKMK4"
    solution = pylie.solve(heavy_top, y0, t_start, t_end, step_length, manifold, method)

    # Verify that the solution is indeed on the manifold
    expexted_norm = np.linalg.norm(y0[3:])
    solution_norm = [np.linalg.norm(solution[3:, i]) for i in range(len(solution.T))]
    for val in solution_norm:
        assert_almost_equal(val, expexted_norm)
    print("Passed test, plotting ...")

    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")
    ax.plot(solution[0, :], solution[1, :], solution[2, :])
    plt.show()

Available numerical schemes

• "E1": Explicit Euler, 1st order
• "RKMK4": Runge-Kutta Munthe-Kaas 4, 4th order