Supplementary information (comments on demo implementations) for Continuation and bifurcation in Nonlinear PDEs – Algorithms, applications, and experiments

General comments. In [6] we consider three example problems for numerical continuation and bifurcation analysis in nonlinear PDEs. Here we comment on the associated pde2path implementations, in the demos seco, sh35disk, and schnakdc, and additionally on the Allen–Cahn dead core demo acdc, which prepares schnakdc, all included in the download DMVdemos at [9]. We start with general comments on pde2path, see also [7, Chapter 5], and then give overviews and specifics for the individual demos. In Table 1 we list acronyms used in [6] and here.

with u = u(x,t) ∈ R N (N components), x ∈ Ω ⊂ R d some bounded domain, d ∈ {1, 2, 3} (the 1D, 2D and 3D case, respectively), and time t ≥ 0, and where (1) can be complemented with various BCs. We refer to [6] or [7] for the meaning of terms in G(u, λ ). The basic idea is to first discretize (1) in space using some FEM built into pde2path, and then treat the resulting high-dimensional ODEs, respectively algebraic systems in case of steady states of (1). As usual in pde2path we assume that all problem data is contained in the MATLAB struct p. This includes the object p.pdeo (with sub-objects fem and grid), which provides methods to generate FEM meshes and assemble FEM matrices M (mass matrix) and K (e.g., Laplacian, including the BCs), and several more, for instance for mesh adaptation. Typical initializations and first continuation steps run as follows (Listings 1-3 contain concrete examples): 1. Call p=stanparam() to initialize the fields in p with defaults values. 2. Call a pdeo constructor, for instance p.pdeo=stanpdeo1D(p,aux) (which discretizes an interval), where here and in the following aux or . . . stands for variable arguments. 3. Initialize p.u with a first solution (or a solution guess, to be corrected in a Newton loop), and append the parameters at the end of p.u. 4. In a function oosetfemops (in the working directory 1 ), use p.pdeo.assema to generate p.mat.M 1 MATLAB always searches the working directory first, which is an easy way to overload pde2path library functions (the dynamical mass matrix, always required) and p.mat.K (a Laplacian), and possibly further FEM matrices, e.g., for BCs. 5. Use p.mat.M and p.mat.K in a function r=sG(p,u) to encode the (discretized) PDE (and the Jacobian in Gu=sGjac(p,u)). Here, the input argument u contains the "PDE unknowns" u and the parameters (appended at the end), and p is typically useful for simple coding, in particular the subfields of p.mat such as the preassembled matrices M and K. 6. Call p=cont(p) to (attempt to) continue the initial solution in some parameter, including bifurcation detection, localization, saving to disk, and plotting. 7. Call p=swibra(dir,bpt,newdir) (or p0=qswibra(dir,bpt,aux); p=seltau(...), see [ to plot bifurcation diagrams, and plotsol(dir,pt,aux) to plot sample solutions. Steps 1-3, and a call to oosetfemops are typically combined into an init-function, for instance p=secoinit(...) in Listing 2.
Remark 1. a) The right hand side G, Jacobian ∂ u G, and a number of further functions needed/used to run pde2path on a problem p, are interfaced via function handles in p.fuha. For instance, you can give the function encoding G in (1) any name, e.g., myrhs, with signature r=myrhs(p,u), and then set p.fuha.sG=@myrhs. In most demos, we simply keep the "standard names" sG and sGjac and code these functions in the respective demo directory. For other handles in p.fuha there are standard choices which we seldomly modify, e.g. p.fuha.savefu=@stansavefu. Functions for which the "default choice" is more likely to be modified include, e.g., • p.fuha.outfu=@stanbra; % signature out=stanbra(p,u). Here the user defines what data (typically parameter values and some norm(s) of u) is used for branch output, e.g., for later plotting of BDs. • p.fuha.lss=@lss; % signature [x,p]=lss(A,b,p). Linear systems solver The default lss is just an interface to MATLAB's \; other options include, e.g., lssbel (bordered elimination) and lssAMG (preconditioned GMRES using ilupack [2]). b) An important feature of the FEM used to spatially discretize the PDEs is its flexibility with respect to the domain shape and BCs, and its flexibility and well established procedures for adaptive mesh refinement. pde2path comes with a number of convenience functions to discretize standard domains (intervals, rectangles, cuboids, disks, sectors, balls, . . . ), and with methods for mesh adaptation in 1D, 2D, and 3D, based on standard error-estimators. Here, we only use some mesh adaptation in the demo acdc, and refer to [7, Chapters 4 and 6] for details.
The demo seco. In [6, §3.1] we consider the 2-component reaction system set up in [4] as a model for semiconductors. Throughout we fix (τ, d) = (8, 0.05), and initially also α = 0.02, and take j 0 as the primary continuation/bifurcation parameter. For all ( j 0 , α), (2) has the unique spatially homogeneous steady state For suitable parameter choices, there are codimension-2 points, near which u * may undergo either a (steady) Turing or a Hopf bifurcation, see [6,Fig. 4].  Fig. 5 and Fig. 6]. secoinit initialization of problem struct p with standard parameter values, call of stanpdeo1D to generate a 1D PDE object (interval, with mesh), initialization of u with the homogeneous steady state, call of oosetfemops to generate the FEM matrices. oosetfemops assemble and store the mass matrix M, and the (1-component) Neumann-Laplacian K. sG,sGjac rhs of (4), and Jacobian. nodaljac "local" derivatives (terms in Jac without spatial derivatives), called in sGjac and spufu. bpjac implements ∂ u (G u φ ) for BP continuation, see also hpjac for HP continuation. secobra mod of library function hobra; subtract steady state from solution for branch output getss convenience function to compute steady state u * from parameters. spufu "spectral" user function, used to plot dispersion relations. Table 2 list the files used to implement and run (2). in the demo seco, and Listing 1 shows the three basic functions essentially needed in a typical pde2path demo. In secoinit we set up the domain, the initial solution (3), based on the parameters passed to secoinit in the script file cmds1, and set a few pde2path control parameters as appropriate for this problem. This, to some extent, is trial and error; moreover, switches and numerical controls such as p.nc.dsmax can be (and often are) changed later any time. Additionally, in line 5 we switch on the bordered elimination linear system solver lssbel, which in 1D, due to the band structure of M and K, usually yields a significant speedup. Here the border of the Jacobian ∂ u H of the extended system H = (G, p), with p the arclength equation [6, (6)], only has width 1 from p, and the 0 in setbel means that ∂ u G itself has no border.
Listing 3 gives the first 9 lines of the script file cmds1, of altogether 70 lines, which contain many plotting commands and somewhat detailed comments. The file is organized in cells, started by %%, which can be executed individually, which we strongly recommend. In the second cell, we use initeig to compute a guess for shifts for Hopf eigenvalues, cf.  Fig.5]. We choose (initial) parameters, the domain size 2l x = 2n w π/k c with n w = 8 and critical wave number k c such that 8 waves fit into Ω at criticality, and choose 40 points per wave for the spatial discretization. Then we compute (via initeig) a guess for the temporal wave number ω 1 near which we aim to detect Hopf eigenvalues, and continue the trivial branch in j 0 . Subsequently, we switch to a number of Turing modes, to the localized Turing mode, and from the generated snaking branch to some localized Hopf modes. Then we deal with plotting.
The further files from Table 2 are optional and/or for special tasks. The function bpjac.m implements ∂ u (G u φ ) which appears in BP continuation, cf. [6, §2.4], and to hpjac.m implementing a similar (large) matrix of derivatives of G u φ r and G u φ i needed for HP continuation, cf.[7, §3.6.1 and §4.4]. In 1D, these functions can be omitted using numerical Jacobians, but in 2D (or even 3D) this becomes rather slow. Finally, secobra is a modification of the library function hobra, which generates [pars; m1; m2; m3; m4] for branch-output, where m 1 = T (period, 0 for steady states), m 2 = max(u 1 ), m 3 = min(u 1 ), and m 4 = u 2 as defined in [6, (51)].
Listing 4 shows the first step to compute the Turing-resp. Hopf bifurcation lines in j 0 -α parameter space ([6, Fig. 4]) via BP and HP continuation. Here we choose α as the new primary active parameter (at position 2 in the parameter vector), and hence j 0 becomes a secondary active parameter. Our main comment is that to improve robustness we here need to increase the parameter p.nc.del (steplength for the finite differences for parameter derivatives, which are always done numerically) from its default 10 −4 to 10 −2 . The HP continuation works similarly via p=hpcontini('0','hpt1',2,'hpc1'), and then we plot the results. Listing 5 first shows a typical call of hoswibra and subsequent continuation of the periodic orbit (first 4 lines). Then we illustrate the rather special task for the system (2), of how to splice together a Turing mode and a Hopf mode to obtain a "mixed solution", see [6, Fig.6].
The demo sh35disk. In [6, §3.1], following [11], we consider the cubic-quintic Swift-Hohenberg (SH) equation ∂ t u = −(1 + ∆) 2 u + εu + νu 3 − u 5 on a disk with Neumann BCs for u and ∆u. Setting u 1 = u and u 2 = ∆u we rewrite the 4th order SH equation as the 2nd order system with a singular dynamical mass matrix M d = 1 0 0 0 , f (u 1 ) = νu 3 1 − u 5 1 , and Neumann BCs for u 1 and u 2 . The implementation again follows standard principles of pde2path demos, see Table 3 for the used files and comments. Compared to the demo seco, the functions q*.m are new. They implementphase conditions [6, §2.4]: After standard preparations in shinit and oosetfemops, qy and qyder are the y-phase conditions for localized solutions on the half disk, see Listing 6, while on the full disk we also need qx and qrot and combinations thereof. Basic script, generating [6,. shinit initialization; setup of half-disk domain and mesh, initialization with u ≡ 0. oosetfemops assembly of mass matrix M, Laplacian K, and of matrices for phase conditions. sG, sGjac rhs of (4), and Jacobian. qy, qyder phase condition q y (u) = ∂ y u old , u = (Ky * p.u) * u, with ∂ u q y (u) = (Ky * p.u) . Similar for ∂ x -or ∂ φ -phase conditions (qx or qrot), where q rot (u) = (−y∂ x + x∂ y )u old , u = (Krot * p.u) * u, and combinations such as qxy for x-and y-phase condition, qxyr (all three) and their derivatives. qyon convenience function to switch on the y-phase condition; similar for qxon and qroton etc. spjac implements ∂ u (G u φ ) for fold continuation. shJ, shbra shJ computes the energy F [6, (53)], and is called and put onto output branch in shbra (modification of library function stanbra). h2fdisk (mirror u from) half disk to full disk. See also h2fdisk0 where instead of mirroring we extend u by 0. Listing 6: Vertical shift phase condition qy and derivative qyder.

Remark 2.
A special feature, also compared to the examples from [7], is that for (4) we use a piecewise quadratic FEM following [5]. The main advantage is that this is more robust than the default P 1 -FEM wrt "branch jumping". See also [8], and [10] for an online version of this demo.
The demo acdc. To show how dead core (DC) problems of type [6, (56) and (63)] can be treated with pde2path, we first consider modifications of the Allen-Cahn type equation from [1], given by where f (u, x) = u(1 − u)(u − a(x)) with a ∈ C 1 ([0, 1], (0, 1)). For fixed x, the ODE d dt u = f (x, u) has the two stable fixed points u = 0 and u = 1, and the unstable fixed point u = a(x). The idea to allow an x-dependent a(x) is that for non-constant choices of a, e.g., periodic a, and small c 2 > 0 there exist stable states u for (5) for which u(x) ≈ 0 on some intervals I 1 , . . . , I n , and u(x) ≈ 1 roughly on the complements, with narrow interfaces in between. This is partly motivated by applications, see [1], but mostly by the fact [3] that for any constant a ≡ a 0 ∈ (0, 1) states with interfaces are unstable and the only stable steady states are the constant states u ≡ 0 or u ≡ 1. Similar results also hold for the case of Dirichlet BCs: For instance, requiring u(0,t) = 1, u(1,t) = 0, the only stable steady state is monotonously decreasing and thus has exactly one interface.
In Fig. 1(a) we show a BD and some sample solutions for the case Ω = (0, 1), (γ, c 2 ) = (0.85, 0.01). Depending on the parameters, in particular γ and c 2 , also some higher bifurcations (with kernels φ = sin( πx), ≥ 2) from the trivial branch lead to branches which develop stable interfaces and subsequently stable dead cores consisting of several disjoint intervals.  (6), Jacobian, and the x-(resp. r-) dependent function a. nloop, nloopext modifications of the Newton loop algorithms, incorporating (9). mypsol customized plotting to visualize DCs, see also mpsol2D. Table 4 list the files for the implementation of (6), in 1D and 2D. The two main issues we need to deal with are the non-differentiability of f from (7) at u = 0, and with keeping u ≥ 0 throughout, in particular during Newton iterations such as and similar for the extended systems [6, (22b)]. To keep u ≥ 0 we make local copies of the pde2path library functions nloop and nloopext in the working directory, and, e.g., in nloop.m add the command u1(1 : p.nu) = max(u1(1 : p.nu), 0) after computing the update u1=u n+1 , and similar in nloopext.m. For u 0 we have ∂ u f (u, λ , x) ∞ and f is not differentiable at u = 0. Thus, to approximate the (non-existing) "Jacobian" for with a small δ > 0, namely δ = 10 −6 . Together, (9) and (10) can be seen as an ad hoc Newton method for (7), which exploits that DCs are valid solutions, and since we do not modify (7), we solve the original problem. 4 In Fig. 2 we essentially study the same problem as in Fig. 1, now in 2D. We use an almost disk shaped domain, which however we slightly disturb to an ellipse with major axis e = 1.1 to break the rotational invariance. We also use a larger diffusion c 2 = 0.05 to avoid too steep interfaces, and start with a rather coarse mesh of nt ≈ 3300 elements. Naturally, steep interfaces suggest adaptive mesh refinement. In pde2path, mesh adaptation is based on user defined elements-to-refine-selector functions, linked as p.fuha.e2rs. The default choice is the library function e2rs, using a standard element wise error estimate η τ , see [7, §4.2.1]. Here, instead of η τ we just use the (discrete) curvature (|∆u|) τ as an ad hoc selector, yielding the refined meshes in Fig. 2(b), with around 10000 elements. This mesh adaptation is needed to continue the nontrivial branches to larger λ , i.e., λ > 20, say. The primary bifurcating branch (c1) shows an elliptic DC for λ > λ c ≈ 9.4 (with free boundary roughly parallel to the domain boundary ∂ Ω), and the second branch shows a DC for λ > 21, which becomes non-convex for λ > 22.5. The demo schnakdc. In [6, §3.3] we consider the Schnakenberg reaction diffusion system ∂ t u = ∆u − u 1/m The term −u 1/m + replaces the standard term −u, and for m > 1 is singular as u 0. This leads to DC pattern formation, including Hopf bifurcations from solution branches with DCs, where the (time-)oscillations at least near bifurcation are strongly localized at the "live" part u > 0. For the implementation we closely follow the demo acdc, where additionally we also use modified Newton loops for periodic orbits, namely honloopexp, similar as in (9). Table 5 lists and comments on the files used in the demo schnakdc.  (11), and Jacobian. nloop, nloopext modifications of the Newton loops, incorporating upos, and similar in honloopext. mypsol customized plotting to visualize DCs, see also mpsol2D, and myhopl (for Hopf orbits).