# pacopy
[](https://pypi.org/project/pacopy/)
[](https://pypi.org/project/pacopy/)
[](https://github.com/nschloe/pacopy)
[](https://pypistats.org/packages/pacopy)
[](https://discord.gg/hnTJ5MRX2Y)
pacopy provides various algorithms of [numerical parameter
continuation](https://en.wikipedia.org/wiki/Numerical_continuation) for ODEs and PDEs in
Python.
pacopy is backend-agnostic, so it doesn't matter if your problem is formulated with
[NumPy](https://numpy.org/), [SciPy](https://scipy.org/),
[FEniCS](https://fenicsproject.org/), [pyfvm](https://github.com/nschloe/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.
Install with
```
pip install pacopy
```
To get started, take a look at the examples below.
Some pacopy documentation is available
[here](https://pacopy.readthedocs.io/en/latest/?badge=latest).
### Examples
#### Basic scalar example
<img src="https://nschloe.github.io/pacopy/simple.svg" width="30%">
Let's start off with a problem where the solution space is scalar. We try to
solve `sin(x) - lambda` for different values of `lambda`, starting at 0.
```python
import math
import matplotlib.pyplot as plt
import pacopy
class SimpleScalarProblem:
def inner(self, a, b):
"""The inner product in the problem domain. For scalars, this is just a
multiplication.
"""
return a * b
def norm2_r(self, a):
"""The norm in the range space; used to determine if a solution has been found."""
return a**2
def f(self, u, lmbda):
"""The evaluation of the function to be solved"""
return math.sin(u) - lmbda
def df_dlmbda(self, u, lmbda):
"""The function's derivative with respect to the parameter. Used in Euler-Newton
continuation.
"""
return -1.0
def jacobian_solver(self, u, lmbda, rhs):
"""A solver for the Jacobian problem. For scalars, this is just a division."""
return rhs / math.cos(u)
problem = SimpleScalarProblem()
lmbda_list = []
values_list = []
def callback(k, lmbda, sol):
# Use the callback for plotting, writing data to files etc.
lmbda_list.append(lmbda)
values_list.append(sol)
pacopy.euler_newton(
problem,
u0=0.0,
lmbda0=0.0,
callback=callback,
max_steps=20,
newton_tol=1.0e-10,
verbose=False,
)
# plot solution
plt.plot(values_list, lmbda_list, ".-")
plt.xlabel("$u_1$")
plt.ylabel("$\\lambda$")
plt.show()
```
#### Simple 2D problem
<img src="https://nschloe.github.io/pacopy/simple2d.svg" width="30%">
A similarly simple example with two unknowns and a parameter. The inner product and
Jacobian solver are getting more interesting.
```python
import matplotlib.pyplot as plt
import numpy as np
import pacopy
import scipy
from mpl_toolkits.mplot3d import Axes3D
class SimpleProblem2D:
def inner(self, a, b):
return np.dot(a, b)
def norm2_r(self, a):
return np.dot(a, a)
def f(self, u, lmbda):
return np.array(
[
np.sin(u[0]) - lmbda - u[1] ** 2,
np.cos(u[1]) * u[1] - lmbda,
]
)
def df_dlmbda(self, u, lmbda):
return -np.ones_like(u)
def jacobian_solver(self, u, lmbda, rhs):
A = np.array(
[
[np.cos(u[0]), -2 * u[1]],
[0.0, np.cos(u[1]) - np.sin(u[1]) * u[1]],
]
)
return np.linalg.solve(A, rhs)
problem = SimpleProblem2D()
# Initial guess and initial parameter value
u0 = np.zeros(2)
lmbda0 = 0.0
# Init lists
lmbda_list = []
values_list = []
def callback(k, lmbda, sol):
lmbda_list.append(lmbda)
values_list.append(sol)
pacopy.euler_newton(problem, u0, lmbda0, callback, max_steps=50, newton_tol=1.0e-10)
# plot solution
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(*np.array(values_list).T, lmbda_list, ".-")
ax.set_xlabel("$u_1$")
ax.set_ylabel("$u_2$")
ax.set_zlabel("$\\lambda$")
# plt.show()
plt.savefig("simple2d.svg", transparent=True, bbox_inches="tight")
plt.close()
```
#### Bratu
<img src="https://nschloe.github.io/pacopy/bratu1d.png" width="30%">
Let's deal with an actual PDE, the classical [Bratu
problem](https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation)
in 1D with Dirichlet boundary conditions. Now, the solution space isn't scalar, but a
vector of length `n` (the values of the solution at given points on the 1D interval).
We now need numpy and scipy, the inner product and Jacobian solver are more complicated.
```python
import matplotlib.pyplot as plt
import numpy as np
import pacopy
import scipy.sparse
import scipy.sparse.linalg
# This is the classical finite-difference approximation
class Bratu1d:
def __init__(self, n):
self.n = n
h = 1.0 / (self.n - 1)
self.H = np.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
)
def inner(self, a, b):
"""The inner product of the problem. Can be np.dot(a, b), but factoring in
the mesh width stays true to the PDE.
"""
return np.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 np.dot(a, a)
def f(self, u, lmbda):
"""The evaluation of the function to be solved"""
out = self.A.dot(u) - lmbda * np.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 = -np.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 * np.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)
problem = Bratu1d(51)
# Initial guess and parameter value
u0 = np.zeros(problem.n)
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)
# use the norm of the currentsolution for plotting on the y-axis
values_list.append(np.sqrt(problem.inner(sol, sol)))
ax1.plot(lmbda_list, values_list, "-x", color="C0")
ax1.set_xlim(0.0, 4.0)
ax1.set_ylim(0.0, 6.0)
plt.close()
# Natural parameter continuation
# pacopy.natural(problem, u0, lmbda0, callback, max_steps=100)
pacopy.euler_newton(
problem, u0, lmbda0, callback, max_steps=500, max_num_retries=10, newton_tol=1.0e-10
)
```
#### Ginzburg–Landau
https://user-images.githubusercontent.com/181628/146639709-90b6e6aa-48ba-418d-9aa4-ec5754f95b93.mp4
The [Ginzburg-Landau
equations](https://en.wikipedia.org/wiki/Ginzburg%E2%80%93Landau_theory) model
the behavior of extreme type-II superconductors under a magnetic field. The
above example (to be found in full detail
[here](tests/test_ginzburg_landau.py)) shows parameter continuation in the
strength of the magnetic field. The plot on the right-hand side shows the
complex-valued solution using [cplot](https://github.com/nschloe/cplot).
Raw data
{
"_id": null,
"home_page": "",
"name": "pacopy",
"maintainer": "",
"docs_url": null,
"requires_python": ">=3.8",
"maintainer_email": "",
"keywords": "mathematics,numerical continuation,numerical methods,differential equations,pde",
"author": "",
"author_email": "Nico Schl\u00f6mer <nico.schloemer@gmail.com>",
"download_url": "",
"platform": null,
"description": "# pacopy\n\n[](https://pypi.org/project/pacopy/)\n[](https://pypi.org/project/pacopy/)\n[](https://github.com/nschloe/pacopy)\n[](https://pypistats.org/packages/pacopy)\n\n[](https://discord.gg/hnTJ5MRX2Y)\n\npacopy provides various algorithms of [numerical parameter\ncontinuation](https://en.wikipedia.org/wiki/Numerical_continuation) for ODEs and PDEs in\nPython.\n\npacopy is backend-agnostic, so it doesn't matter if your problem is formulated with\n[NumPy](https://numpy.org/), [SciPy](https://scipy.org/),\n[FEniCS](https://fenicsproject.org/), [pyfvm](https://github.com/nschloe/pyfvm), or any\nother Python package. The only thing the user must provide is a class with some simple\nmethods, e.g., a function evaluation `f(u, lmbda)`, a Jacobian a solver\n`jacobian_solver(u, lmbda, rhs)` etc.\n\nInstall with\n\n```\npip install pacopy\n```\n\nTo get started, take a look at the examples below.\n\nSome pacopy documentation is available\n[here](https://pacopy.readthedocs.io/en/latest/?badge=latest).\n\n### Examples\n\n#### Basic scalar example\n\n<img src=\"https://nschloe.github.io/pacopy/simple.svg\" width=\"30%\">\n\nLet's start off with a problem where the solution space is scalar. We try to\nsolve `sin(x) - lambda` for different values of `lambda`, starting at 0.\n\n```python\nimport math\nimport matplotlib.pyplot as plt\nimport pacopy\n\n\nclass SimpleScalarProblem:\n def inner(self, a, b):\n \"\"\"The inner product in the problem domain. For scalars, this is just a\n multiplication.\n \"\"\"\n return a * b\n\n def norm2_r(self, a):\n \"\"\"The norm in the range space; used to determine if a solution has been found.\"\"\"\n return a**2\n\n def f(self, u, lmbda):\n \"\"\"The evaluation of the function to be solved\"\"\"\n return math.sin(u) - lmbda\n\n def df_dlmbda(self, u, lmbda):\n \"\"\"The function's derivative with respect to the parameter. Used in Euler-Newton\n continuation.\n \"\"\"\n return -1.0\n\n def jacobian_solver(self, u, lmbda, rhs):\n \"\"\"A solver for the Jacobian problem. For scalars, this is just a division.\"\"\"\n return rhs / math.cos(u)\n\n\nproblem = SimpleScalarProblem()\n\nlmbda_list = []\nvalues_list = []\n\n\ndef callback(k, lmbda, sol):\n # Use the callback for plotting, writing data to files etc.\n lmbda_list.append(lmbda)\n values_list.append(sol)\n\n\npacopy.euler_newton(\n problem,\n u0=0.0,\n lmbda0=0.0,\n callback=callback,\n max_steps=20,\n newton_tol=1.0e-10,\n verbose=False,\n)\n\n# plot solution\nplt.plot(values_list, lmbda_list, \".-\")\nplt.xlabel(\"$u_1$\")\nplt.ylabel(\"$\\\\lambda$\")\nplt.show()\n```\n\n#### Simple 2D problem\n\n<img src=\"https://nschloe.github.io/pacopy/simple2d.svg\" width=\"30%\">\n\nA similarly simple example with two unknowns and a parameter. The inner product and\nJacobian solver are getting more interesting.\n\n```python\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pacopy\nimport scipy\nfrom mpl_toolkits.mplot3d import Axes3D\n\n\nclass SimpleProblem2D:\n def inner(self, a, b):\n return np.dot(a, b)\n\n def norm2_r(self, a):\n return np.dot(a, a)\n\n def f(self, u, lmbda):\n return np.array(\n [\n np.sin(u[0]) - lmbda - u[1] ** 2,\n np.cos(u[1]) * u[1] - lmbda,\n ]\n )\n\n def df_dlmbda(self, u, lmbda):\n return -np.ones_like(u)\n\n def jacobian_solver(self, u, lmbda, rhs):\n A = np.array(\n [\n [np.cos(u[0]), -2 * u[1]],\n [0.0, np.cos(u[1]) - np.sin(u[1]) * u[1]],\n ]\n )\n return np.linalg.solve(A, rhs)\n\n\nproblem = SimpleProblem2D()\n# Initial guess and initial parameter value\nu0 = np.zeros(2)\nlmbda0 = 0.0\n\n# Init lists\nlmbda_list = []\nvalues_list = []\n\n\ndef callback(k, lmbda, sol):\n lmbda_list.append(lmbda)\n values_list.append(sol)\n\n\npacopy.euler_newton(problem, u0, lmbda0, callback, max_steps=50, newton_tol=1.0e-10)\n\n# plot solution\nfig = plt.figure()\nax = fig.add_subplot(111, projection=\"3d\")\nax.plot(*np.array(values_list).T, lmbda_list, \".-\")\nax.set_xlabel(\"$u_1$\")\nax.set_ylabel(\"$u_2$\")\nax.set_zlabel(\"$\\\\lambda$\")\n# plt.show()\nplt.savefig(\"simple2d.svg\", transparent=True, bbox_inches=\"tight\")\nplt.close()\n```\n\n#### Bratu\n\n<img src=\"https://nschloe.github.io/pacopy/bratu1d.png\" width=\"30%\">\n\nLet's deal with an actual PDE, the classical [Bratu\nproblem](https://en.wikipedia.org/wiki/Liouville%E2%80%93Bratu%E2%80%93Gelfand_equation)\nin 1D with Dirichlet boundary conditions. Now, the solution space isn't scalar, but a\nvector of length `n` (the values of the solution at given points on the 1D interval).\nWe now need numpy and scipy, the inner product and Jacobian solver are more complicated.\n\n```python\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pacopy\nimport scipy.sparse\nimport scipy.sparse.linalg\n\n\n# This is the classical finite-difference approximation\nclass Bratu1d:\n def __init__(self, n):\n self.n = n\n h = 1.0 / (self.n - 1)\n\n self.H = np.full(self.n, h)\n self.H[0] = h / 2\n self.H[-1] = h / 2\n\n self.A = (\n scipy.sparse.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(self.n, self.n))\n / h**2\n )\n\n def inner(self, a, b):\n \"\"\"The inner product of the problem. Can be np.dot(a, b), but factoring in\n the mesh width stays true to the PDE.\n \"\"\"\n return np.dot(a, self.H * b)\n\n def norm2_r(self, a):\n \"\"\"The norm in the range space; used to determine if a solution has been found.\"\"\"\n return np.dot(a, a)\n\n def f(self, u, lmbda):\n \"\"\"The evaluation of the function to be solved\"\"\"\n out = self.A.dot(u) - lmbda * np.exp(u)\n out[0] = u[0]\n out[-1] = u[-1]\n return out\n\n def df_dlmbda(self, u, lmbda):\n \"\"\"The function's derivative with respect to the parameter. Used in Euler-Newton\n continuation.\n \"\"\"\n out = -np.exp(u)\n out[0] = 0.0\n out[-1] = 0.0\n return out\n\n def jacobian_solver(self, u, lmbda, rhs):\n \"\"\"A solver for the Jacobian problem.\"\"\"\n M = self.A.copy()\n d = M.diagonal().copy()\n d -= lmbda * np.exp(u)\n M.setdiag(d)\n # Dirichlet conditions\n M.data[0][self.n - 2] = 0.0\n M.data[1][0] = 1.0\n M.data[1][self.n - 1] = 1.0\n M.data[2][1] = 0.0\n return scipy.sparse.linalg.spsolve(M.tocsr(), rhs)\n\n\nproblem = Bratu1d(51)\n# Initial guess and parameter value\nu0 = np.zeros(problem.n)\nlmbda0 = 0.0\n\nlmbda_list = []\nvalues_list = []\n\n\ndef callback(k, lmbda, sol):\n # Use the callback for plotting, writing data to files etc.\n fig = plt.figure()\n ax1 = fig.add_subplot(111)\n ax1.set_xlabel(\"$\\\\lambda$\")\n ax1.set_ylabel(\"$||u||_2$\")\n ax1.grid()\n\n lmbda_list.append(lmbda)\n # use the norm of the currentsolution for plotting on the y-axis\n values_list.append(np.sqrt(problem.inner(sol, sol)))\n\n ax1.plot(lmbda_list, values_list, \"-x\", color=\"C0\")\n ax1.set_xlim(0.0, 4.0)\n ax1.set_ylim(0.0, 6.0)\n plt.close()\n\n\n# Natural parameter continuation\n# pacopy.natural(problem, u0, lmbda0, callback, max_steps=100)\n\npacopy.euler_newton(\n problem, u0, lmbda0, callback, max_steps=500, max_num_retries=10, newton_tol=1.0e-10\n)\n```\n\n#### Ginzburg\u2013Landau\n\nhttps://user-images.githubusercontent.com/181628/146639709-90b6e6aa-48ba-418d-9aa4-ec5754f95b93.mp4\n\nThe [Ginzburg-Landau\nequations](https://en.wikipedia.org/wiki/Ginzburg%E2%80%93Landau_theory) model\nthe behavior of extreme type-II superconductors under a magnetic field. The\nabove example (to be found in full detail\n[here](tests/test_ginzburg_landau.py)) shows parameter continuation in the\nstrength of the magnetic field. The plot on the right-hand side shows the\ncomplex-valued solution using [cplot](https://github.com/nschloe/cplot).\n",
"bugtrack_url": null,
"license": "",
"summary": "Numerical continuation in Python",
"version": "0.2.8",
"project_urls": {
"Homepage": "https://github.com/sigma-py/pacopy",
"Issues": "https://github.com/sigma-py/pacopy/issues"
},
"split_keywords": [
"mathematics",
"numerical continuation",
"numerical methods",
"differential equations",
"pde"
],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "cbe175ad649c19705108d41b99f35800d2586f31322dfd93147cb7e0770a5598",
"md5": "9f2a838920dd4ee326fb089e8407a55f",
"sha256": "8baccc6b64ffb1946f68c910cece959563fed5667fe7b120132b8f5e1bdc96ca"
},
"downloads": -1,
"filename": "pacopy-0.2.8-py3-none-any.whl",
"has_sig": false,
"md5_digest": "9f2a838920dd4ee326fb089e8407a55f",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": ">=3.8",
"size": 22710,
"upload_time": "2024-03-15T08:33:06",
"upload_time_iso_8601": "2024-03-15T08:33:06.445609Z",
"url": "https://files.pythonhosted.org/packages/cb/e1/75ad649c19705108d41b99f35800d2586f31322dfd93147cb7e0770a5598/pacopy-0.2.8-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2024-03-15 08:33:06",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "sigma-py",
"github_project": "pacopy",
"travis_ci": false,
"coveralls": false,
"github_actions": false,
"lcname": "pacopy"
}
JSON
