SciPy

Fundamental algorithms for scientific computing in Python.

Link	Web site
Language	Python
License	BSD-3-Clause license
Documentation	documentation
Git repository	@scipy

Lorenz equations

We would like to solve the Lorenz equations, a classical chaotic system given by:

\[\begin{split} \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - x z \\ \dot{z} &= x y - \beta z \end{cases} \end{split}\]

with parameter \(\sigma=10\), \(\rho = 28\) and \(\beta = \frac{8}{3}\), and the initial state \((x_0, y_0, z_0) = (1,1,1)\). We solve it with a classical Runge-Kutta method of order 4.

To solve Lorenz equations with a classical 4th order Runge-Kutta RK(4,4) we need to define it:

class RK44(scipy.integrate._ivp.rk.RungeKutta):
    order = 4
    n_stages = 4
    A = np.array([
        [0, 0, 0, 0],
        [1/2, 0, 0, 0],
        [0, 1/2, 0, 0],
        [0, 0, 1, 0]
    ])
    B = np.array([1/6, 1/3, 1/3, 1/6])
    C = np.array([0, 1/2, 1/2, 1])

and now call it with scipy.integrate.solve_ivp function:

sol = solve_ivp(lorenz_system, t_span, y0, method=RK44, first_step=dt)

To solve a problem with SciPy we first write our equation as a ODE of the form:

\[ \dot{u} = f(t, u) \]

and the user provide the function \(f\) as:

  def f(t, u):
    # ...
    return du

where u the current state of the function, t the current time and output du \(f(t,u)\).

def lorenz_system(t, y):
    sigma = 10.
    rho = 28.
    beta = 8./3.

    return np.asarray([
        sigma * (y[1] - y[0]),
        y[0]*(rho - y[2]) - y[1],
        y[0] * y[1] - beta * y[2]
    ])

After chose a method, we can solve the problem between initial time and final time

sol = solve_ivp(lorenz_system, t_span, y0, method=RK44,
                first_step=dt, max_step=dt, rtol=1.0, atol=1.0)

this function returns an object with solution at each time (and also dense output properties).

For the complet example, see lorenz.py source file.

Transport equation

In this example we would like to solve the following PDE:

\[ \partial_t u + a \partial_x u = 0 \]

with \(t>0\), on the torus \(x\in[0, 1)\), a velocity \(a=1\) and the initial condition given by a hat function:

\[\begin{split} u(0, x) = \begin{cases} x - 0.25 & \text{if } x\in[0.25, 0.5[ \\ -x + 0.75 & \text{if } x\in[0.5, 0.75[ \\ 0 & \text{else} \end{cases} \end{split}\]

We choose a first order up-wind scheme to estimate the \(x\) derivative and a forward Euler method for the time discretization.

SciPy provides mainly adaptive time step methods, so there is no classical 4th order Runge-Kutta RK(4,4) nor explicit Euler (or forward Euler or RK(1, 1)), we need to define it:

class Euler(scipy.integrate._ivp.rk.RungeKutta):
    order = 1
    n_stages = 1
    A = np.array([
        [0]
    ])
    B = np.array([1])
    C = np.array([0])

and now call it with scipy.integrate.solve_ivp function:

sol = solve_ivp(upwind, t_span, y0, method=Euler, first_step=dt)

We define the up-wind scheme as:

def upwind(t, y):
    dy = np.zeros_like(y)

    dy[0] = - (np.max([a, 0]) * (y[0] - y[-1]) +
               np.min([a, 0]) * (y[1] - y[0])) / dx
    dy[1:-1] = -(np.max([a, 0]) * (y[1:-1] - y[:-2]) +
                 np.min([a, 0]) * (y[2:] - y[1:-1])) / dx
    dy[-1] = - (np.max([a, 0]) * (y[-1] - y[-2]) +
                np.min([a, 0]) * (y[0] - y[-1])) / dx

    return dy

The time loop is the same as for Lorenz equation.

For the complet example, see transport.py source file.

Arenstorf orbit

The Arenstorf orbit problem is a classical problem to test adaptive time step methods:

\[\begin{split} \begin{cases} \ddot{x} &= x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} &= y - 2\dot{x} - \frac{}{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases} \end{split}\]

with initial condition \((x,\dot{x},y,\dot{y})=(0.994, 0, 0, -2.001585106)\), \(r_1\) and \(r_2\) given by

\[ r_1 = \sqrt{(x+\mu)^2 + y^2},\quad r_2 = \sqrt{(x-1+\mu)^2 + y^2} \]

and with parameter \(\mu = 0.012277471\).

First of all, we need to rewrite this problem into a first order derivative equation in time

\[\begin{split} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix}, \qquad \begin{cases} \dot{y}_1 = y_3 \\ \dot{y}_2 = y_4 \\ \dot{y}_3 = y_1 + 2y_4 - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ \dot{y}_4 = y_2 - 2y_3 - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \\ \end{cases} \end{split}\]

We define this system as:

def arenstorf_system(t, y):
    mu = 0.012277471

    y1 = y[0]
    y2 = y[1]
    y3 = y[2]
    y4 = y[3]

    r1 = np.sqrt((y1 + mu) * (y1 + mu) + y2 * y2)
    r2 = np.sqrt((y1 - 1. + mu) * (y1 - 1. + mu) + y2 * y2)

    dy1 = y3
    dy2 = y4
    dy3 = y1 + 2 * y4 - (1 - mu) * (y1 + mu) / \
        (r1 * r1 * r1) - mu * (y1 - 1 + mu) / (r2 * r2 * r2)
    dy4 = y2 - 2 * y3 - (1 - mu) * y2 / (r1 * r1 * r1) - \
        mu * y2 / (r2 * r2 * r2)

    return np.asarray([dy1, dy2, dy3, dy4])

We solve this example with given method RK45 which is the method RK5(4) 7M in [DP80] (mainly call DOPRI5) and DOP853 which is the method RK8(7) 13M in [PD81] (mainly call DOPRI8).

The time loop is the same as for Lorenz equation, for RK45 method

sol = solve_ivp(arenstorf_system, t_span, y0, method="RK45",
                first_step=dt, rtol=1e-5, atol=1e-5)

and for DOP853 method

sol = solve_ivp(arenstorf_system, t_span, y0, method="DOP853",
                first_step=dt, rtol=1e-5, atol=1e-5)

For the complet example, see arenstorf.py source file.