bice package

Subpackages

Module contents

BICE: Bifurcation Continuation Engine.

A numerical path continuation and bifurcation analysis package written in Python. It provides tools for tracking branches of solutions to nonlinear algebraic equations, identifying bifurcation points, and performing time-stepping for evolutionary PDEs.

class BifurcationDiagram

Bases: object

Manages a collection of branches.

Includes methods for plotting.

Initialize the BifurcationDiagram.

active_branch: Branch: storage for the currently active branch

branches: list[Branch]: list of branches

current_solution() → Solution

Return the latest solution in the diagram.

Returns:: The latest solution of the active branch.
Return type:: Solution

get_branch_by_ID(branch_id: int) → Branch | None

Return a branch by its ID.

Parameters:: branch_id – The ID of the branch.
Returns:: The branch object or None if not found.
Return type:: Branch or None

load_branch(filename: str) → None

Load a branch from a file into the diagram.

The file should have been stored with Branch.save(filename).

Parameters:: filename – The filename to load from.

new_branch(active: bool = True) → Branch

Create a new branch.

Parameters:: active – Whether to set the new branch as active.
Returns:: The newly created branch.
Return type:: Branch

norm_name: name of the norm

parameter_name: str | None: name of the continuation parameter

plot(ax: Any) → None

Plot the bifurcation diagram.

Parameters:: ax – The matplotlib axes to plot onto.

remove_branch_by_ID(branch_id: int) → None

Remove a branch from the BifurcationDiagram by its ID.

Parameters:: branch_id – The ID of the branch to remove.

xlim: tuple[float, float] | None: x-limits of the diagram

ylim: tuple[float, float] | None: y-limits of the diagram

class Branch

Bases: object

Stores a list of solution objects from a parameter continuation.

Initialize the Branch.

add_solution_point(solution: Solution) → None

Add a solution to the branch.

Parameters:: solution – The solution object to add.

bifurcations() → list[Solution]

List all bifurcation points on the branch.

Returns:: List of solutions that are bifurcation points.
Return type:: list

data(only: str | None = None) → tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]

Return the list of parameters and norms of the branch.

Parameters:: only – Optional filter: “stable”, “unstable”, or “bifurcations”.
Returns:: Tuple of (parameter values, norm values).
Return type:: tuple[RealArray, RealArray]

id: unique identifier of the branch

is_empty() → bool: Check if the current branch is empty.

norm_vals() → ndarray[tuple[Any, ...], dtype[float64]]

List of solution norm values along the branch.

Returns:: Array of norm values.
Return type:: RealArray

parameter_vals() → ndarray[tuple[Any, ...], dtype[float64]]

List of continuation parameter values along the branch.

Returns:: Array of parameter values.
Return type:: RealArray

remove_solution_point(solution: Solution) → None

Remove a solution from the branch.

Parameters:: solution – The solution object to remove.

save(filename: str) → None

Store the branch to the disk.

Parameters:: filename – The filename to save to.

solutions: list[Solution]: list of solutions along the branch

class EigenSolver

Bases: object

A wrapper to the powerful iterative Eigensolver ARPACK.

Finds eigenvalues and eigenvectors of an eigenproblem.

Initialize the EigenSolver.

latest_eigenvalues: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]] | None: results of the latest eigenvalue computation

latest_eigenvectors: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]] | None: results of the latest eigenvector computation

shift: The shift used for the shift-invert method in the iterative eigensolver. If shift != None, the eigensolver will find the eigenvalues near the value of the shift first

Solve the eigenproblem A*x = v*x for the eigenvalues v and the eigenvectors x.

If an unsigned integer k is given, the iterative eigensolver ARPACK will calculate the k first eigenvalues sorted by the largest real part. Otherwise a direct eigensolver will be used to calculate all eigenvalues.

If a mass matrix M is given, the generalized eigenvalue A*x = v*M*x will be solved.

Parameters:

A – The system matrix.
M – The mass matrix (optional).
k – Number of eigenvalues to compute (optional).

Returns:

A tuple containing (eigenvalues, eigenvectors).

Return type:

tuple[Array, Array]

tol: convergence tolerance of the eigensolver

class Equation(shape: int | tuple[int, ...] | None = None)

Bases: object

Base class for algebraic (Cauchy) equations.

Holds equations of the form M du/dt = rhs(u, t, r) where M is the mass matrix, u is the vector of unknowns, t is the time and r is a parameter vector. This may include ODEs and PDEs. All custom equations must inherit from this class and implement the rhs(u) method. Time and parameters are implemented as member attributes.

This is a very fundamental class. Specializations of the Equation class exist for covering more intricate types of equations, i.e., particular discretizations for spatial fields, e.g., finite difference schemes or pseudospectral methods.

An equation has a ‘shape’ that should at all times be equal to the shape of the unknowns ‘u’.

Initialize the equation.

Parameters:: shape – The shape of the unknowns. If None, defaults to ().

adapt() → tuple[float, float] | None

Adapt the equation to the solution (mesh refinement or similar).

May be overridden for specific types of equations, do not forget to adapt Equation.u_history as well!

group: EquationGroup: Reference to a group of equations that this equation belongs to. Initialized with the DummyEquationGroup that only contains the current equation.

is_coupled: Does the equation couple to any other unknowns? If it is coupled, then all unknowns and methods of this equation will have the full dimension of the problem and need to be mapped to the equation’s variables accordingly. Otherwise, they only have the dimension of this equation.

jacobian(u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → ndarray | spmatrix

Calculate the Jacobian J = d rhs(u) / du for the unknowns u.

Defaults to automatic calculation of the Jacobian using finite differences.

Parameters:: u – The vector of unknowns.
Returns:: The Jacobian matrix.
Return type:: Matrix

load(data: dict[str, Any]) → None

Restore unknowns / parameters / etc. from the given dictionary.

The dictionary was created by Equation.save(). Equation.load() is the inverse of Equation.save(). May be overridden for loading more stuff for specific types of equations.

Parameters:: data – The data dictionary to load from.

mass_matrix() → ndarray | spmatrix

Return the mass matrix M.

The mass matrix determines the linear relation of the rhs to the temporal derivatives: M * du/dt = rhs(u).

Returns:: The mass matrix.
Return type:: Matrix

property ndofs: int

Return the total number of unknowns / degrees of freedom of the equation.

Returns:: The total number of degrees of freedom.
Return type:: int

plot(ax: Any) → None

Plot the solution into a matplotlib axes object.

Parameters:: ax – The matplotlib axes to plot into.

reshape(shape: int | tuple[int, ...]) → None

Change the shape of the equation / the equation’s unknowns.

Parameters:: shape – The new shape of the unknowns.

rhs(u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Calculate the right-hand side of the equation 0 = rhs(u).

Parameters:: u – The vector of unknowns.
Returns:: The result of the right-hand side evaluation.
Return type:: Array
Raises:: NotImplementedError – If the method is not implemented by the subclass.

save() → dict[str, Any]

Save everything that is relevant for this equation to a dict.

The Problem class will call this and save the dict to the disk. May be overridden for saving more stuff for specific types of equations.

Returns:: A dictionary containing the state of the equation.
Return type:: dict

property shape: tuple[int, ...]

self.u.shape.

Returns:: The shape of the unknowns.
Return type:: tuple[int, …]
Type:: Return the shape of the equation’s unknowns

u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]: The equation’s storage for the unknowns

u_history: list[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]]: a history of the unknowns, needed e.g. for implicit schemes

class EquationGroup(equations: list[Equation | EquationGroup] | None = None, auto_add: bool = True)

Bases: object

Groups multiple equations into a single new equation (a system of equations).

All properties and functions are assembled from the subequations. EquationGroups may even form hierarchical trees, where one group of equations serves as a subequation to another one.

Initialize the EquationGroup.

Parameters:

equations – A list of equations or equation groups to include.
auto_add – Whether to automatically add the given equations.

__repr__() → str

Return a string representation of the EquationGroups in the terminal.

Returns:: The string representation.
Return type:: str

add_equation(eq: Equation | EquationGroup) → None

Add an equation to the group.

Parameters:: eq – The equation or equation group to add.

equations: list[Equation | EquationGroup]: the list of sub-equations (or even sub-groups-of-equations)

group: EquationGroup | None: optional reference to a parent EquationGroup

idx: dict[Equation | EquationGroup, slice]: The indices of the equation’s unknowns to the group’s unknowns and vice versa

property is_coupled: bool

Return whether the group is coupled.

A group of equations should never couple to other groups.

Returns:: Always False.
Return type:: bool

jacobian(u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → ndarray | spmatrix

Calculate the Jacobian J = d rhs(u) / du for the unknowns u.

Parameters:: u – The vector of unknowns.
Returns:: The Jacobian matrix.
Return type:: Matrix

list_equations() → list[Equation]

Return a flattened list of all equations in the group and sub-groups.

Returns:: The list of equations.
Return type:: list[Equation]

map_unknowns() → None

Create the mapping from equation unknowns to group unknowns.

This ensures that group.u[idx[eq]] = eq.u.ravel() where idx is the mapping.

mass_matrix() → ndarray | spmatrix

Return the mass matrix.

The mass matrix determines the linear relation of the rhs to the temporal derivatives: M * du/dt = rhs(u).

Returns:: The mass matrix.
Return type:: Matrix

property ndofs: int

Return the number of unknowns / degrees of freedom of the group.

Returns:: The total number of degrees of freedom.
Return type:: int

remove_equation(eq: Equation | EquationGroup) → None

Remove an equation from the group.

Parameters:: eq – The equation or equation group to remove.

rhs(u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Calculate the right-hand side of the group 0 = rhs(u).

Parameters:: u – The vector of unknowns.
Returns:: The result of the right-hand side evaluation.
Return type:: Array

property shape: tuple[int, ...]

Return the shape of the unknowns.

Returns:: The shape of the unknowns vector.
Return type:: tuple[int, …]

property u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

combined unknowns of the sub-equations.

Returns:: The combined vector of unknowns.
Return type:: Array
Type:: Return the unknowns of the system

class MyNewtonSolver

Bases: AbstractNewtonSolver

Reference implementation of a simple ‘text book’ Newton solver.

Initialize the AbstractNewtonSolver.

solve(f: Callable[[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], u0: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]], jac: Callable[[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], ndarray | spmatrix]) → ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Solve the system using a custom Newton iteration.

Parameters:

f – The function to find the root of.
u0 – The initial guess.
jac – The Jacobian of the function.

Return type:

The solution vector.

class NewtonSolver

Bases: AbstractNewtonSolver

A Newton solver that uses scipy.optimize.root for solving.

The method (algorithm) to be used can be adjusted with the attribute ‘method’. NOTE: does not work with sparse Jacobians, but converts it to a dense matrix instead!

Initialize the NewtonSolver.

method: choose from the different methods of scipy.optimize.root NOTE: method = “krylov” might be faster, but then we can use NewtonKrylovSolver directly

solve(f: Callable[[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], u0: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]], jac: Callable[[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]], ndarray | spmatrix] | None = None) → ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Solve using scipy.optimize.root.

Parameters:

f – The function.
u0 – Initial guess.
jac – The Jacobian (optional).

Return type:

The solution vector.

class Problem

Bases: object

Base class for all algebraic problems.

It is an aggregate of (one or many) governing algebraic equations, initial and boundary conditions, constraints, etc. Also, it provides all the basic properties, methods and routines for treating the problem, e.g., time-stepping, solvers or plain analysis of the solution. Custom problems should be implemented as children of this class.

Initialize the Problem.

adapt() → None: Adapt the problem/equations to the solution (e.g. by mesh refinement).

add_equation(eq: Equation | EquationGroup) → None

Add an equation to the problem.

Parameters:: eq – The equation to add.

bifurcation_diagram: The bifurcation diagram of the problem holds all branches and their solutions

continuation_parameter: tuple[object, str] | None: The continuation parameter is defined by passing an object and the name of the object’s attribute that corresponds to the continuation parameter as a tuple

continuation_step() → None: Perform a parameter continuation step.

continuation_stepper: The continuation stepper for parameter continuation

eigen_solver: The eigensolver for eigenvalues and -vectors

eq: None | Equation | EquationGroup: the equation (or system of equation) that governs the problem

generate_bifurcation_diagram(parameter_lims: tuple[float, float] = (-1000000000.0, 1000000000.0), norm_lims: tuple[float, float] = (-1000000000.0, 1000000000.0), max_recursion: int = 4, max_steps: float = 1000000000.0, detect_circular_branches: bool = True, ax: Any | None = None, plotevery: int = 30) → None

Automatically generate a full bifurcation diagram within the given bounds.

Branch switching will be performed automatically up to the given maximum recursion level.

Parameters:

parameter_lims – Limits for the continuation parameter (min, max).
norm_lims – Limits for the norm (min, max).
max_recursion – Maximum recursion depth for sub-branch continuation.
max_steps – Maximum number of steps per branch.
detect_circular_branches – Stop when solution reaches starting point again.
ax – Axes object to live plot the diagram.
plotevery – Plotting frequency (every N steps).

get_continuation_parameter() → float

Return the value of the continuation parameter.

Returns:: The current value of the parameter.
Return type:: float

history: The history of the unknown values is accessed and managed with the Problem.history object

jacobian(u: ArrayLike) → ndarray | spmatrix

Calculate the Jacobian of the system J = d rhs(u) / du for the unknowns u.

Parameters:: u – The vector of unknowns.
Returns:: The Jacobian matrix.
Return type:: Matrix

list_equations() → list[Equation]

List all equations that are part of the problem.

Returns:: The list of equations.
Return type:: list[Equation]

load(data: dict[str, Any] | Solution | str) → None

Load the current solution from the given data.

Problem.load(…) is the inverse of Problem.save(…).

Parameters:: data – Filename, Solution object, or dictionary.

locate_bifurcation(ev_index: int | None = None, tolerance: float = 1e-05) → bool

Locate the closest bifurcation using bisection method.

Finds point where the real part of the eigenvalue closest to zero vanishes.

Parameters:

ev_index – Optional index of the eigenvalue that corresponds to the bifurcation.
tolerance – Threshold at which the value is considered zero.

Returns:

True if the location converged, False otherwise.

Return type:

bool

locate_bifurcation_using_constraint(eigenvector: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → None

Locate the bifurcation of the given eigenvector.

Parameters:: eigenvector – The eigenvector to use for the constraint.

log(*args: Any, **kwargs: Any) → None

Wrap print() for log messages.

Log messages are printed only if verbosity is switched on.

Parameters:

*args – Positional arguments for print.
**kwargs – Keyword arguments for print.

mass_matrix() → ndarray | spmatrix

Return the mass matrix.

The mass matrix determines the linear relation of the rhs to the temporal derivatives: M * du/dt = rhs(u).

Returns:: The mass matrix.
Return type:: Matrix

property ndofs: int

Return the number of unknowns / degrees of freedom of the problem.

Returns:: The number of degrees of freedom.
Return type:: int

new_branch() → None: Create a new branch in the bifurcation diagram and prepare for continuation.

newton_solve() → None

Solve the system rhs(u) = 0 for u with Newton’s method.

Updates self.u with the solution.

newton_solver: The Newton solver for finding roots of equations

norm() → float

Return the default norm of the solution.

Used for bifurcation diagrams. Defaults to the L2-norm of the unknowns.

Returns:: The norm value.
Return type:: float

Plot everything to the given axes.

Axes may be given explicitly or as a list of axes, that is then expanded. The plot may include the solution of the equations, the bifurcation diagram, the eigenvalues and the eigenvectors.

Parameters:

sol_ax – Axes for plotting the solution.
bifdiag_ax – Axes for plotting the bifurcation diagram.
eigvec_ax – Axes for plotting the eigenvectors.
eigval_ax – Axes for plotting the eigenvalues.

remove_equation(eq: Equation | EquationGroup) → None

Remove an equation from the problem.

Parameters:: eq – The equation to remove.

rhs(u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]) → ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Calculate the right-hand side of the system 0 = rhs(u).

Parameters:: u – The vector of unknowns.
Returns:: The residuals vector.
Return type:: Array

save(filename: str | None = None) → dict[str, Any]

Save the current solution to the file <filename>.

Returns a dictionary of the serialized data.

Parameters:: filename – Optional filename to save to.
Returns:: The serialized data.
Return type:: dict

set_continuation_parameter(val: float | floating) → None

Set the value of the continuation parameter.

Parameters:: val – The new value for the parameter.

settings: The settings (tolerances, switches, etc.) are held by this ProblemSettings object

solve_eigenproblem() → tuple[ndarray[tuple[Any, ...], dtype[float64 | complexfloating]], ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]]

Calculate the eigenvalues and eigenvectors of the Jacobian.

The method will only calculate as many eigenvalues as requested with self.settings.neigs.

Returns:: Tuple of (eigenvalues, eigenvectors).
Return type:: tuple[Array, Array]

switch_branch(ev_index: int | None = None, amplitude: float = 0.001, locate: bool = True) → bool

Attempt to switch branches in a bifurcation.

Parameters:

ev_index – The index of the eigenvalue corresponding to the bifurcation.
amplitude – The amplitude of the perturbation.
locate – Whether to attempt to locate the bifurcation point first.

Returns:

True if successful, False otherwise.

Return type:

bool

time: float: Time variable

time_step() → None: Integrate in time with the assigned time-stepper.

time_stepper: TimeStepper: The time-stepper for integration in time

property u: ndarray[tuple[Any, ...], dtype[float64 | complexfloating]]

Getter for unknowns of the problem.

Returns:: The flattened vector of unknowns.
Return type:: Array

class Profiler

Bases: object

Static class for accessing/controlling the profiling of the code.

execution_times: dict[str, float] = {}

static is_active() → bool: Check if the Profiler is active/running.

static print_summary(nested: bool = True) → None

Print a summary on the execution times of the decorated methods.

Parameters:: nested – Whether to show a nested tree view or a flattened list.

static start() → None: (Re)start the Profiler.

class Solution(problem: Problem | None = None)

Bases: object

Stores the solution of a problem.

Includes the relevant parameters and some information on the solution.

Initialize the Solution.

Parameters:: problem – The problem instance to store the state from.

bifurcation_type(update: bool = False) → str

Return the type of bifurcation.

Parameters:: update – Force update of the cached type.
Returns:: The bifurcation type string (e.g., “BP”, “HP”, “+”, “-“).
Return type:: str

branch: Branch | None: optional reference to the corresponding branch

data: dict[str, Any]: The current problem state as a dictionary of data (equation’s unknowns and parameters)

get_neighboring_solution(distance: int) → Solution | None

Get access to the neighboring solution in the branch.

Parameters:: distance – The index distance (can be negative).
Returns:: The neighboring solution or None if not found/out of bounds.
Return type:: Solution or None

id: unique identifier of the solution

is_bifurcation() → bool

Check if the solution point is a bifurcation.

Returns:: True if it is a bifurcation point.
Return type:: bool

is_stable() → bool | None

Check if the solution is stable.

Returns:: True if stable, False if unstable, None if unknown.
Return type:: bool or None

property neigenvalues_crossed: int | None

Return how many eigenvalues have crossed the imaginary axis.

Returns:: The number of eigenvalues crossed, or None if unknown.
Return type:: int or None

property nimaginary_eigenvalues_crossed: int | None

Return how many eigenvalues have crossed the imaginary axis.

Returns:: The number of eigenvalues crossed, or None if unknown.
Return type:: int or None

norm: float: value of the solution norm

nunstable_eigenvalues: int | None: number of true positive eigenvalues

nunstable_imaginary_eigenvalues: int | None: number of true positive and imaginary eigenvalues

p: float

profile(method: Callable) → Callable

Decorate a function to trigger profiling of its execution time.

Parameters:: method – The function to profile.
Returns:: The wrapped function.
Return type:: Callable