PaperlayerView exact revision on arXiv

Source-linked AI summary

pde2path - A Matlab package for continuation and bifurcation in 2D elliptic systems

Hannes Uecker, Daniel Wetzel, Jens D. M. Rademacher

arXiv:1208.3112v2math.APmath.NA

TL;DR

pde2path addresses the lack of accessible general-purpose continuation and bifurcation tools for two-dimensional elliptic PDE systems. It combines Matlab finite-element problem definitions with one-parameter continuation, bifurcation handling, examples, and customization. The examples demonstrate continuation and bifurcation computations, while the package remains limited to one parameter and simple bifurcations.

  • Problem

    Few general continuation and bifurcation tools work out of the box for non-expert users studying two-dimensional elliptic PDE systems.

  • Method

    pde2path uses a modular Matlab finite-element framework in which users define PDE coefficients, Jacobians, geometry, boundary conditions, and initial guesses, then run continuation and bifurcation workflows.

  • Results

    The examples include detailed continuation and bifurcation computations for systems such as Schnakenberg reaction–diffusion and Rayleigh–Bénard convection.

  • Takeaways & Limitations

    pde2path provides an out-of-the-box starting point for exploring the rich solution and bifurcation structures of two-dimensional elliptic systems.

  • Takeaways & Limitations

    The current package supports only one-parameter continuation and simple bifurcations, with undesired branch switching also occurring in some Schnakenberg computations.

Abstract

from arXiv · show

pde2path is a free and easy to use Matlab continuation/bifurcation package for elliptic systems of PDEs with arbitrary many components, on general two dimensional domains, and with rather general boundary conditions. The package is based on the FEM of the Matlab pdetoolbox, and is explained by a number of examples, including Bratu's problem, the Schnakenberg model, Rayleigh-Benard convection, and von Karman plate equations. These serve as templates to study new problems, for which the user has to provide, via Matlab function files, a description of the geometry, the boundary conditions, the coefficients of the PDE, and a rough initial guess of a solution. The basic algorithm is a one parameter arclength continuation with optional bifurcation detection and branch-switching. Stability calculations, error control and mesh-handling, and some elementary time-integration for the associated parabolic problem are also supported. The continuation, branch-switching, plotting etc are performed via Matlab command-line function calls guided by the AUTO style. The software can be downloaded from www.staff.uni-oldenburg.de/hannes.uecker/pde2path, where also an online documentation of the software is provided such that in this paper we focus more on the mathematics and the example systems.

1 Introduction

pde2path addresses the scarcity of general, accessible continuation and bifurcation tools for two-dimensional elliptic PDE systems. It provides a flexible Matlab-based framework with command-line workflows, examples, customization, and several analysis capabilities.

  • pde2path targets elliptic PDE systems in two spatial dimensions, where few general continuation and bifurcation tools work out of the box for non-experts.
  • The package represents systems as G(u, λ) := −∇· (c ⊗∇u) + au −b ⊗∇u −f = 0 on bounded two-dimensional domains.
  • Its supported generalized Neumann boundary conditions include zero flux and stiff-spring approximations of Dirichlet conditions.
  • The package also supports parabolic time integration for generating continuation initial conditions and Jacobian eigenvalue calculations for stability inspection.
  • Users specify geometry, boundary conditions, PDE coefficients, and a rough initial guess, then call Matlab functions for continuation, bifurcation detection, branch switching, and time integration.
  • The implementation emphasizes modular, transparent data structures so users can customize problems and add features, favoring readability over speed.
  • Current limitations include Matlab and pdetoolbox requirements, slower performance than some Fortran FEM implementations, one-parameter continuation, and detection of only simple-eigenvalue bifurcations.

2 Some basics of continuation and bifurcation

The paper presents pseudo-arclength continuation as an extended system combining the PDE residual with a tangent-based hyperplane condition. Weighting, Newton correction, eigenvalue monitoring, and branch-switching support continuation through folds and simple bifurcations, with documented scope limits.

  • 2.1 Arclength continuation: Pseudo-arclength continuation augments G(u, λ) = 0 with a tangent-based scalar condition, allowing solution branches to be followed using predictors and Newton correctors.
  • 2.1 Arclength continuation: The method can use a chord corrector that avoids repeated Jacobian evaluation at the cost of a usually modest increase in iterations.
  • 2.1 Arclength continuation: The weight ξ controls the relative emphasis on changes in u versus λ and must be adapted when parameter scaling changes.
  • 2.1 Arclength continuation: Each continuation step predicts along the normalized tangent, corrects by Newton iteration with adaptive step size, and recomputes the tangent.
  • 2.1 Arclength continuation: Bifurcation detection uses sign changes of detA, implicitly targeting simple eigenvalue crossings, while nearby eigenvalues support stability and consistency checks.
  • 2.1 Arclength continuation: The one-parameter version does not specially treat folds, excludes bifurcations involving even numbers of eigenvalues, and detects only simple bifurcations.
  • 2.1 Arclength continuation: After locating a bifurcation by bisection, branch switching constructs a tangent for the alternate branch and resumes continuation with chosen ξ and step size.
  • 2.2 Switching back and forth to the natural parametrization: The natural-parameterized Newton scheme is usually slightly faster than arclength continuation, but arclength typically locates bifurcation points better on nearly horizontal branches.

3 Some scalar problems in pde2path

The tutorial introduces pde2path through basic examples, with directory-based demos that provide commands and verbose scripts for running and adapting problems.

  • The tutorial begins with basic examples whose bracketed names identify corresponding subdirectories in the demos directory.
  • Each demo includes command scripts containing example commands and comments intended as templates for typical pde2path usage.
  • The examples are organized to support learning the package by running and adapting existing problem setups.

3.1 Bratu’s problem (bratu)

The Bratu example shows how to formulate, initialize, continue, and analyze a two-dimensional elliptic problem in pde2path. It also illustrates Jacobian choices, mesh considerations, bifurcation structure, and solver behavior.

  • Problem setup: The Bratu problem is a scalar elliptic equation, −∆u − f(u, λ) = 0 with f(u, λ) = −10(u − λe^u), on the unit square with zero-flux boundary conditions.
  • Analytical structure: Its homogeneous solution branch starts at (0, 0), while simple bifurcation points are predicted from the Laplacian eigenvalues and can lead to further bifurcations on secondary branches.
  • Representation and setup: A problem is represented by a structure p containing coefficient and boundary-condition function handles, geometry and mesh fields, the current solution and tangent, and algorithm controls.
  • Representation and setup: Users define coefficient and Jacobian functions, initialize geometry, mesh, parameters, and a starting point, then run continuation and branch switching through pde2path functions.
  • Mesh handling: Delaunay-based and rectangular-domain meshes are available, but the authors recommend experimenting with both because Delaunay meshes can provide more robust numerics after refinement.
  • Jacobians and numerical accuracy: Analytically assembled Jacobians are faster by orders of magnitude on large meshes, whereas numerical Jacobians can improve Newton convergence but may mask poor FEM approximations on underresolved meshes.
  • Linear solvers: For scalar problems, iterative solvers can outperform direct solvers when np > 10^5, while the supplied iterative routines are intended mainly as templates for problem-specific solvers.

3.2 The Allen–Cahn equation with Dirichlet boundary conditions (ac)

The section demonstrates pde2path on a cubic–quintic Allen–Cahn equation with Dirichlet conditions, using continuation in λ and subsequent continuation in μ. It also introduces mesh adaptation and semi-implicit time integration for perturbed solutions.

  • Continuation in λ: Choosing Lx = 1, Ly = 0.9 and μ = 0.25 breaks square symmetry and produces simple bifurcation points while continuing in λ.The first bifurcation parameters are λ11 = 1.3784, λ21 = 3.2289, and λ12 = 3.6630.
  • Numerical handling: After branch switching, mesh adaptation is activated every five continuation steps to support the computed branches.The section directs readers to acdemo.m and accmds.m for continuation, solution perturbation, and time-integration examples.
  • Bifurcation structure: Figure 4 reports an elementary bifurcation diagram with no secondary bifurcations and branch mode structures fixed at bifurcation.Panels (b) and (c) show selected branch points, while panel (d) shows a solution obtained by continuing μ from one selected point.
  • Time integration: pde2path provides semi-implicit Euler time integration by assembling M + hK(u^(n)) and solving for u^(n+1) at each time step.Perturbing an unstable solution can lead subsequent trajectories toward stable solutions on the unstable manifold example.
  • Time integration: The current time integrator uses an elementary linear-system solve at every step rather than special-case matrix reuse or factorization.The implementation notes that optimization is possible when coefficients are independent of u, but is not currently used.

3.3 The Allen–Cahn equation with mixed λ-dependent boundary conditions (achex)

This example extends the Allen–Cahn problem to a hexagonal domain with parameter-dependent mixed boundary conditions. pde2path performs continuation, bifurcation detection, branch switching, and error estimation in this setting.

  • Problem setup: The Allen–Cahn equation is reconsidered on a hexagonal domain with mixed Dirichlet–Neumann boundary conditions.The domain is based on a square whose top boundary is shifted, with ΓD and ΓN partitioning the boundary.
  • Boundary conditions: The boundary conditions impose n · ∇u = 0 on ΓN and u = λx on ΓD.The Dirichlet condition is parameter dependent, while the Neumann condition is homogeneous.
  • Implementation: The implementation uses a stiff-spring approximation through gnbcs to impose the Dirichlet boundary condition.The supplied commands set a penalty scale of 10^4 for the Dirichlet part and zero data for the Neumann part.
  • Numerical results: Starting from the trivial zero solution, continuation produces the bifurcation diagram in Figure 5, with bifurcation detection and branch switching functioning in the example.The reported error estimate remains below 0.01.
  • Numerical results: Both directions of the r branch are generated by branch switching, continuation, loading the resulting point, reversing the step direction, and continuing again.The two continuation runs use opposite ds values to proceed in the other direction.

3.4 A quasilinear Allen–Cahn equation (acql)

The quasilinear Allen–Cahn example modifies the diffusion coefficient to depend on u, requiring a more complicated Jacobian. This change breaks the equation’s sign symmetry and changes the first bifurcation type.

  • Problem setup: The quasilinear equation replaces constant diffusion with the coefficient 0.25 + δu + γu^2 in the divergence-form operator.It is posed with homogeneous Dirichlet conditions on the same rectangular domain.
  • Jacobian: The linearization contains diffusion, reaction, Laplacian, and gradient terms depending on u, δ, and γ.The resulting operator is given explicitly as Gu(u,λ)v in the section.
  • Jacobian: The implementation incorporates the additional derivatives into fu, b111, and b112 using mesh-derived gradients and Laplacians.The coefficient construction uses pdegrad and pdeprtni in acqljac.m.
  • Bifurcation changes: The δu diffusion term breaks the u 7→−u symmetry while leaving the bifurcation points from the trivial branch unchanged.The bifurcations themselves change despite the unchanged locations.
  • Bifurcation changes: The first bifurcation changes from pitchfork to transcritical when the symmetry-breaking diffusion term is introduced.Figure 6 illustrates the resulting branch structure and solution plots for δ = −0.2 and γ = 0.05.

3.5 An Allen–Cahn equation with global coupling (acgc)

This example treats an Allen–Cahn equation with global feedback, which introduces a nonlocal term beyond the local PDE framework. The FEM Jacobian is represented as a rank-one modification and solved using a Sherman–Morrison formula.

  • Problem setup: The global-coupling problem adds the feedback term −λ⟨u⟩ to a fixed-parameter Allen–Cahn equation on a square domain.Here μ = 0.1 and the equation uses homogeneous Dirichlet conditions.
  • Nonlocal structure: Global coupling is nonlocal because ⟨u⟩ integrates u over Ω, so it does not fit the framework assuming a local f.The text relates such coupling to surface catalysis, semiconductors, gas discharges, and shadow systems.
  • FEM Jacobian: Ignoring the global-coupling term in continuation works only for small |λ| and can destroy Newton convergence for larger |λ|.The implementation therefore retains the nonlocal contribution in the linear solve.
  • FEM Jacobian: At the FEM level, the Jacobian becomes Gu(u) = K − λνη^T, where ν and η encode the discretized global-average operation.The vector η is formed from triangle areas and interpolation from nodal values to triangle values.
  • Linear solver: The Sherman–Morrison formula solves the rank-one-modified system without explicitly forming the full matrix K − λνη^T.The resulting expression uses solves involving K^-1 and a scalar correction factor α.
  • Implementation: The implementation uses customized solvers and mesh-dependent global variables ν and η, recalculating them when mesh refinement requires it.The relevant routines include acgcf.m, acgcjac.m, gclss, and gcblss.
  • Solution structure: Positive global feedback raises the plateau of selected solutions above the zero of f(u) = u + u^3 − u^5, while negative feedback can push boundary values below zero.Other branches contain somewhat localized solutions, and the problem is symmetric under (u,λ) 7→(−u,λ).

3.6 First summary, and some remarks on customization

pde2path emphasizes a modular, command-line workflow that supports customization while retaining accessible defaults and standard continuation operations.

  • Customization: Users can maintain multiple function versions and customize any pde2path function by copying and modifying it for a particular problem.The package limits user-definable functions to a selected set while keeping broader customization possible.
  • Plotting: Plotting supports flexible argument/value lists, with plotbra documented as the main exception to the package’s generally small input/output interfaces.Users are directed to plotbra.m and demo calls for detailed plotting options.
  • Code organization: The software keeps problem-specific calculations separate from structure-modifying functions, with few global variables by default.Functions such as cont, swibra, and meshref modify the problem structure p, including its mesh.
  • Usage: Typical usage proceeds through predefined command-line functions, including continuation, bifurcation detection, branch switching, time integration, and mesh handling.Table 8 summarizes the main function calls and pseudocode for p=cont(p).

4 Some prototype Reaction–Diffusion Systems

The prototype reaction–diffusion examples demonstrate pde2path on chemotaxis and Schnakenberg systems across rectangular and complex domains. They also expose practical continuation issues and motivate parallel multi-continuation near bifurcations.

  • Chemotaxis: The chemotaxis example recovers published bifurcation diagrams on rectangular domains with general cross-diffusion and homogeneous Neumann boundary conditions.The model is encoded through coefficient functions, including a nonsymmetric diffusion tensor depending on u and λ.
  • Schnakenberg model: The Schnakenberg model exhibits many bifurcations on small domains, requiring pmcont and findbif to support parallel continuation and locate initial bifurcations from the homogeneous branch.The model is semilinear with a diagonal constant diffusion matrix.
  • Chemotaxis: The chemotaxis calculation can use assembled Jacobians on a 2376-triangle mesh, while numerical-Jacobian runs require considerably finer meshes and adaptive refinement.Finer meshes reduce error estimates but remove the speed advantage of assembled Jacobians; boundary-gradient evaluation contributes to the difficulty.
  • Chemotaxis: For D = 1/4, r = 1.52, and a 1 × 4 domain, the chemotaxis computations identify bifurcation values with reasonable accuracy and agree well with earlier results.The reported exception is a loop near λ = 20.5 on the (1, 1) branch containing two bifurcations.
  • Schnakenberg model: Standard continuation can switch undesirably from the stripe branch to the beans branch, motivating pmcont as a modification that addresses branch switching near bifurcations.Reducing the continuation parameter only postpones the unwanted switching, while pmcont also incorporates parallel computing for speedup.
  • Patterns and domains: Further branch switching produces additional patterns, including snaking between stripes and hexagons, while the animal-shaped-domain example studies bifurcation directions from a trivial branch.The domain uses Neumann boundary conditions and is set up graphically through the package’s geometry and boundary-condition tools.
  • Continuation strategy: For p.resfac=1 and p.mst=1, pmcont is roughly equivalent to cont; with p.mst> 1, it uses parallel computing and can help avoid convergence problems and undesired branch switching.The two versions remain separate because cont is simpler to hack and implements Keller’s basic algorithm.

5 Three classical examples from physics

pde2path is demonstrated on multiphysics elliptic systems, including condensates, convection, and plate buckling, using continuation, bifurcation detection, stability tools, and mesh handling. The examples expose both the package’s flexibility and numerical limitations in complex two-dimensional systems.

  • Overview: The examples cover multicomponent Bose–Einstein condensates, Rayleigh–Bénard convection, and von Kármán plate buckling with boundary-condition implementations suited to systems with many components.The von Kármán system is presented as the largest and most complicated example in component count and boundary-condition implementation.
  • Bose–Einstein condensates: For condensates, pde2path continues localized multipole solutions in λ = ω, tracking deformation from quadrupoles through azimuthons to vortices.The modulation depth p measures soliton-intensity deformation, and a typical quadrupole continuation uses stiff-spring Dirichlet conditions on a square domain.
  • Bose–Einstein condensates: Phase invariance gives condensate solutions continuous families and a zero Jacobian eigenvalue, although the numerical discretization perturbs this eigenvalue sufficiently for continuation.The example uses mesh refinement and obtains an error estimate below 0.01 for the reported continuation.
  • Rayleigh–Bénard convection: In Rayleigh–Bénard convection, pde2path implements the streamfunction system with analytical Jacobians and compares no-slip and stress-free boundary conditions.The continuation parameter is the Rayleigh number R, with Prandtl number σ set to 1.
  • Rayleigh–Bénard convection: No secondary bifurcations appear up to R = 900 for no-slip conditions, whereas stress-free conditions produce an imperfect pitchfork caused by mesh and boundary-condition asymmetries.Increasing mesh resolution brings the diagram closer to a symmetry-breaking pitchfork, but coarse meshes can shift branch locations by order 100 in R.
  • von Kármán plates: For the regularized von Kármán plate system, the tutorial detects a secondary bifurcation between the first two primary branches, while higher bifurcations require mesh refinement.The tutorial mesh has 1250 triangles; refinement improves the error estimate but makes typical continuation steps take a couple of minutes.

6 Discussion

The discussion positions pde2path as an accessible tool for exploring the rich continuation and bifurcation structure of two-dimensional elliptic systems. It remains a first step, with broader bifurcation handling and genuine multiparameter continuation identified as priorities.

  • Scope and contribution: pde2path addresses the numerical and structural challenges of continuation and bifurcation analysis for two-dimensional elliptic systems, including their rich solution structure.The authors present it as a general tool intended to work essentially out of the box for non-expert users.
  • Open limitations: The package does not yet handle bifurcations through nonsimple eigenvalues, which are common in two-dimensional systems with symmetries.This is listed as a priority for more general bifurcation handling.
  • Open limitations: The current framework is limited to basic one-parameter continuation and simple bifurcations because genuine multiparameter continuation has not yet been implemented.Travelling-wave bifurcations are given as an example that generically require a second parameter and an additional phase condition.
Loading 1208.3112v2…