Skip to main content

Numerical continuation in Python

Project description

pacopy

CircleCI codecov Code style: black Documentation Status PyPi Version GitHub stars

pacopy provides various algorithms of numerical parameter continuation for PDEs in Python.

pacopy is backend-agnostic, so it doesn't matter if your problem is formulated with SciPy, FEniCS, pyfvm, or any other Python package. The only thing the user must provide is a class with some simple methods, e.g., a function evaluation f(u, lmbda), a Jacobian a solver jacobian_solver(u, lmbda, rhs) etc.

Some pacopy documentation is available here.

Examples

Bratu

The classical Bratu problem in 1D with Dirichlet boundary conditions. To reproduce the plot, you first have to specify the problem; this is the classical finite-difference approximation:

class Bratu1d(object):
    def __init__(self):
        self.n = 51
        h = 1.0 / (self.n - 1)

        self.H = numpy.full(self.n, h)
        self.H[0] = h / 2
        self.H[-1] = h / 2

        self.A = (
            scipy.sparse.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(self.n, self.n))
            / h ** 2
        )
        return

    def inner(self, a, b):
        """The inner product of the problem. Can be numpy.dot(a, b), but factoring in
        the mesh width stays true to the PDE.
        """
        return numpy.dot(a, self.H * b)

    def norm2_r(self, a):
        """The norm in the range space; used to determine if a solution has been found.
        """
        return numpy.dot(a, a)

    def f(self, u, lmbda):
        """The evaluation of the function to be solved
        """
        out = self.A.dot(u) - lmbda * numpy.exp(u)
        out[0] = u[0]
        out[-1] = u[-1]
        return out

    def df_dlmbda(self, u, lmbda):
        """The function's derivative with respect to the parameter. Used in Euler-Newton
        continuation.
        """
        out = -numpy.exp(u)
        out[0] = 0.0
        out[-1] = 0.0
        return out

    def jacobian_solver(self, u, lmbda, rhs):
        """A solver for the Jacobian problem.
        """
        M = self.A.copy()
        d = M.diagonal().copy()
        d -= lmbda * numpy.exp(u)
        M.setdiag(d)
        # Dirichlet conditions
        M.data[0][self.n - 2] = 0.0
        M.data[1][0] = 1.0
        M.data[1][self.n - 1] = 1.0
        M.data[2][1] = 0.0
        return scipy.sparse.linalg.spsolve(M.tocsr(), rhs)

Then pass the object to any of pacopy's methods, e.g., the Euler-Newton (arclength) continuation:

problem = Bratu1d()
# Initial guess
u0 = numpy.zeros(problem.n)
# Initial parameter value
lmbda0 = 0.0

lmbda_list = []
values_list = []
def callback(k, lmbda, sol):
    # Use the callback for plotting, writing data to files etc.
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.set_xlabel("$\\lambda$")
    ax1.set_ylabel("$||u||_2$")
    ax1.grid()

    lmbda_list.append(lmbda)
    values_list.append(numpy.sqrt(problem.inner(sol, sol)))

    ax1.plot(lmbda_list, values_list, "-x", color="#1f77f4")
    ax1.set_xlim(0.0, 4.0)
    ax1.set_ylim(0.0, 6.0)
    return

# Natural parameter continuation
# pacopy.natural(problem, u0, lmbda0, callback, max_steps=100)

pacopy.euler_newton(
    problem, u0, lmbda0, callback, max_steps=500, newton_tol=1.0e-10
)

Ginzburg–Landau

ginzburg-landau

The Ginzburg-Landau equations model the behavior of extreme type-II superconductors under a magnetic field. The above example (to be found in full detail here) shows parameter continuation in the strength of the magnetic field. The plot on the right-hand side shows the complex-valued solution using cplot.

Installation

pacopy is available from the Python Package Index, so simply type

pip install -U pacopy

to install or upgrade.

Testing

To run the pacopy unit tests, check out this repository and type

pytest

Distribution

To create a new release

  1. bump the __version__ number,

  2. publish to PyPi and GitHub:

    make publish
    

License

pacopy is published under the MIT license.

Project details


Download files

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

Files for pacopy, version 0.1.1
Filename, size File type Python version Upload date Hashes
Filename, size pacopy-0.1.1-py2.py3-none-any.whl (10.0 kB) File type Wheel Python version py2.py3 Upload date Hashes View hashes
Filename, size pacopy-0.1.1.tar.gz (18.5 kB) File type Source Python version None Upload date Hashes View hashes

Supported by

Elastic Elastic Search Pingdom Pingdom Monitoring Google Google BigQuery Sentry Sentry Error logging AWS AWS Cloud computing DataDog DataDog Monitoring Fastly Fastly CDN DigiCert DigiCert EV certificate StatusPage StatusPage Status page