[go: nahoru, domu]

Skip to content

Python implementation of solvers for differential algebraic equation's (DAEs) that should be added to scipy one day.

License

Notifications You must be signed in to change notification settings

JonasBreuling/scipy_dae

Repository files navigation

scipy_dae - solving differential algebraic equations (DAE's) and implicit differential equations (IDE's) in Python

Actions Status Code coverage status badge License: BSD 3 PyPI

Python implementation of solvers for differential algebraic equations (DAE's) and implicit differential equations (IDE's) that should be added to scipy one day. See the GitHub repository.

Currently, two different methods are implemented.

  • Implicit Radau IIA methods of order 2s - 1 with arbitrary number of odd stages.
  • Implicit backward differentiation formula (BDF) of variable order with quasi-constant step-size and stability/ accuracy enhancement using numerical differentiation formula (NDF).

More information about both methods are given in the specific class documentation.

To pique your curiosity

The Kármán vortex street solved by a finite element discretization of the weak form of the incompressible Navier-Stokes equations using FEniCS and the three stage Radau IIA method.

Karman

Basic usage

The Robertson problem of semi-stable chemical reaction is a simple system of differential algebraic equations of index 1. It demonstrates the basic usage of the package.

import numpy as np
import matplotlib.pyplot as plt
from scipy_dae.integrate import solve_dae


def F(t, y, yp):
    """Define implicit system of differential algebraic equations."""
    y1, y2, y3 = y
    y1p, y2p, y3p = yp

    F = np.zeros(3, dtype=y.dtype)
    F[0] = y1p - (-0.04 * y1 + 1e4 * y2 * y3)
    F[1] = y2p - (0.04 * y1 - 1e4 * y2 * y3 - 3e7 * y2**2)
    F[2] = y1 + y2 + y3 - 1 # algebraic equation

    return F


# time span
t0 = 0
t1 = 1e7
t_span = (t0, t1)
t_eval = np.logspace(-6, 7, num=1000)

# initial conditions
y0 = np.array([1, 0, 0], dtype=float)
yp0 = np.array([-0.04, 0.04, 0], dtype=float)

# solver options
method = "Radau"
# method = "BDF" # alternative solver
atol = rtol = 1e-6

# solve DAE system
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval)
t = sol.t
y = sol.y

# visualization
fig, ax = plt.subplots()
ax.set_xlabel("t")
ax.plot(t, y[0], label="y1")
ax.plot(t, y[1] * 1e4, label="y2 * 1e4")
ax.plot(t, y[2], label="y3")
ax.set_xscale("log")
ax.legend()
ax.grid()
plt.show()

Robertson

Advanced usage

More examples are given in the examples directory, which includes

Install

An editable developer mode can be installed via

python -m pip install -e .[dev]

The tests can be started using

python -m pytest --cov