Detection of high codimensional bifurcations in variational PDEs

TL;DR

This paper develops a numerical method to detect high codimensional bifurcations in parameter-dependent PDEs, using bifurcation test equations and the Infinite-dimensional Splitting Lemma, applicable to both continuous and discretized problems.

Contribution

It introduces a unified approach for detecting A-series bifurcations in variational PDEs, simplifying analysis and computation in Banach spaces.

Findings

Successfully detects swallowtail bifurcations in a Bratu-type problem.

Provides a unified framework applicable to PDEs and their discretizations.

Simplifies the analysis of high codimensional bifurcations in variational PDEs.

Abstract

We derive bifurcation test equations for A-series singularities of nonlinear functionals and, based on these equations, we propose a numerical method for detecting high codimensional bifurcations in parameter-dependent PDEs such as parameter-dependent semilinear Poisson equations. As an example, we consider a Bratu-type problem and show how high codimensional bifurcations such as the swallowtail bifurcation can be found numerically. In particular, our original contributions are (1) the use of the Infinite-dimensional Splitting Lemma, (2) the unified and simplified treatment of all A-series bifurcations, (3) the presentation in Banach spaces, i.e. our results apply both to the PDE and its (variational) discretization, (4) further simplifications for parameter-dependent semilinear Poisson equations (both continuous and discrete), and (5) the unified treatment of the continuous problem and its discretisation.

Equations

$$\{\begin{array}{rl}\Delta u+f(u,\lambda )& =0\\ u{|}_{\partial \Omega }& =0\end{array}$$

$${\mathrm{D}}^{2}S(0)(u,v)=⟨Tu,v{⟩}_E\text{for all }u,v\in E.$$

$$\forall u\in U,\forall v\in \mathrm{rg}T:⟨{\nabla }_{y}S(u),v{⟩}_E=\mathrm{D}S(u)v.$$

$$S(\bar{x},\bar{y})=\frac{1}{2}⟨T\bar{y},\bar{y}⟩+r(\bar{x})$$

$$\begin{array}{r}F(0)=0,\mathrm{D}F(0)=0,\\ \forall x\in U_{\ker T}:{\nabla }_{y}S(x,F(x))=0.\end{array}$$

$$r(x)=S(x,F(x)).$$

$${\frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}|}_{x=0}r(x)=0,k=3,\dots ,n,{\frac{{\mathrm{d}}^{n+1}}{\mathrm{d}{x}^{n+1}}|}_{x=0}r(x)\ne 0.$$

$$\left\{j=(j_1,\dots ,j_{n-k+1}):j_l\in \mathbb{N}\cup \{0\},\sum _{l=1}^{n-k+1}{j}_l=k,\sum _{l=1}^{n-k+1}l\cdot j_l=n\right\}.$$

$$\forall \xi \in \{s\alpha |s\in \mathbb{R}{\}}^{{\perp }_E}:$$

$$\sum _{k=1}^{n-2}\sum _{j\in {\mathcal{J}}_k^{n-2}}\frac{(n-2)!}{j!}{S}_0^{(k+1)}\left(\xi \otimes \underset{l=1}{\overset{n-k-1}{⨂}}\frac{{\left({\frac{{\mathrm{d}}^l}{\mathrm{d}{s}^l}|}_{s=0}(s\alpha +F(s))\right)}^{\otimes j_l}}{{(l!)}^{j_l}}\right)=0.$$

$$\sum _{k=1}^n\sum _{j\in {\mathcal{J}}_k^n}\frac{n!}{j!}{S}_0^{(k)}\left(\underset{l=1}{\overset{n-k+1}{⨂}}\frac{{\left({\frac{{\mathrm{d}}^l}{\mathrm{d}{s}^l}|}_{s=0}(s\alpha +F(s))\right)}^{\otimes j_l}}{{(l!)}^{j_l}}\right)=0.$$

$$\ker T={\mathrm{span}}_{\mathbb{R}}\{\alpha \}$$

$$S_0^{(3)}(\alpha ,\alpha ,\alpha )=0.$$

$$S_0^{(2)}(v,\xi )=-S_0^{(3)}(\alpha ,\alpha ,\xi )\forall \xi \in E$$

$$S_0^{(4)}(\alpha ,\alpha ,\alpha ,\alpha )-3S_0^{(2)}(v,v)=0.$$

$$S_0^{(2)}(w,\xi )=-S_0^{(4)}(\alpha ,\alpha ,\alpha ,\xi )-3S_0^{(3)}(\alpha ,v,\xi )\forall \xi \in E$$

$$S_0^{(5)}({\alpha }^{\otimes 5})-15S_0^{(3)}(\alpha ,v,v)+10S_0^{(3)}(\alpha ,\alpha ,w)=0.$$

$$B_n(x_1,\dots ,x_n)=\sum _{k=1}^n\sum _{j\in {\mathcal{J}}_k^n}\frac{n!}{j!}\prod _{l=1}^{n-k+1}{\left(\frac{x_l}{l!}\right)}^{j_l},$$

$$\begin{array}{rl}{B}_0=& 1,\\ B_1(x_1)=& x_1,\\ B_2(x_1,x_2)=& x_1^2+x_2,\\ B_3(x_1,x_2,x_3)=& x_1^3+3x_1{x}_2+x_3,\\ B_4(x_1,x_2,x_3,x_4)=& x_1^4+6x_1^2{x}_2+4x_1{x}_3+3x_2^2+x_4,\\ B_5(x_1,x_2,x_3,x_4,x_5)=& x_1^5+10x_2{x}_1^3+15x_2^2{x}_1+10x_3{x}_1^2+10x_3{x}_2+5x_4{x}_1+x_5.\end{array}$$

$$\exp \left(\sum _{k=1}^{\infty }\frac{x_k}{k!}{y}^k\right)=\sum _{n=0}^{\infty }\frac{1}{n!}{B}_n(x_1,\dots ,x_n)y^n.$$

$$B_4(x_1,x_2,x_3,x_4)=x_1^4+6x_1^2{x}_2+4x_1{x}_3+3x_2^2+x_4$$

$$r^{(n)}(s)=B_n(\alpha +F^{\prime }(s),F^{\prime \prime }(s),\dots ,F^{(n)}(s)).$$

$$r^{(n)}(s)=\sum _{k=1}^n\sum _{j\in {\mathcal{J}}_k^n}\frac{n!}{j!}{S}^{(k)}(s\alpha +F(s))\left(\underset{l=1}{\overset{n-k+1}{⨂}}\frac{{\left(\frac{{\mathrm{d}}^l}{\mathrm{d}{s}^l}(s\alpha +F(s))\right)}^{\otimes j_l}}{{(l!)}^{j_l}}\right),$$

$$r(0)$$

$$r^{\prime }(0)$$

$$r^{\prime \prime }(0)$$

$$r^{\prime \prime \prime }(0)$$

$$r^{\prime \prime \prime \prime }(0)$$

$$r^{\prime \prime \prime \prime \prime }(0)$$

$$+15S_0^{(3)}(\alpha \otimes (F^{\prime \prime }(0){)}^{\otimes 2})+10S_0^{(2)}(F^{\prime \prime }(0)\otimes F^{\prime \prime \prime }(0)).$$

$$r^{(i)}(0)=B_i(\alpha ,F^{\prime \prime }(0),\dots ,F^{(i-2)}(0),c_2,c_1),$$

$$\ker T={\mathrm{span}}_{\mathbb{R}}\{\alpha \}$$

$$\forall s\in I,\forall \xi \in \mathrm{rg}T:S^{(1)}(s\alpha +F(s))(\xi )=0.$$

Full text

††: Nonlinearity

1 Introduction

Nonlinear systems of ordinary (ODEs) or partial differential equations (PDEs) are the basis for modeling many problems in science and technology, including biological pattern formation, viscous fluid flow phenomena, chemical reactions and crystal growth processes. These mathematical models usually depend on a number of parameters and in terms of applications it is crucial to understand the qualitative dependence of the solution on the model parameters.

To explore the solution set of a parameter-dependent system of ordinary or partial differential equations subject to initial and/or boundary conditions, it is helpful to find parameter values and points in the phase space at which bifurcations occur, i.e. at which the structure of the solution set changes. Most prominent is the fold bifurcation in which two branches of solutions merge and die as one parameter is varied (see its bifurcation diagram on the left in Figure 1). A fold bifurcation is a generic phenomenon, i.e. it persists under small perturbations of the parameter-dependent system of differential equations.

Bifurcations can be classified into low and high codimensional bifurcations according to their codimension where the codimension of a bifurcation is defined as the minimal number of parameters in which that bifurcation type occurs. Low codimensional bifurcations such as folds have been studied extensively in the literature. Detecting such points in the parameter-phase portrait, following the solution branches and classifying them is an important task for getting a better understanding of the physical properties of the underlying model [1]. The detection of bifurcation points and the numerical computation of the associated bifurcation diagram is typically based on following solution branches using numerical continuation methods. Numerical continuation has been an active research area, see [2, 3, 4, 5, 6, 7, 8, 9], for instance, and dynamical systems software packages like AUTO [10, 11], CONTENT [12] and MATCONT [13, 14] are available for the detection of bifurcations in ODEs. Applying continuation methods to PDEs, however, leads to the additional challenge of large systems of equations since these methods are typically based on finite element or finite difference discretizations of PDEs, as in the MATLAB package pde2path [15], for instance. For the continuation process, predictor-corrector algorithms are typically used, i.e. nonlinear systems have to be solved in each continuation step using Newton-type methods [16]. Some convergence results are available [17, 18, 19].

If more than one parameter is present in the definition of the underlying mathematical model, more complicated bifurcations may occur. In particular, the presence of many parameters is required for bifurcation points or singular points of high codimension to occur such that they are persistent under perturbations of the dynamical system. High codimensional bifurcations act like organising centres in bifurcation diagrams, i.e. they determine which bifurcations happen close to singular points [20], and occur in many biological and chemical systems [21, 22, 23]. Therefore, it is of interest to find high codimensional singularities. Just as in detecting low codimensional bifurcations, one typically follows a branch of solutions using numerical continuation methods [5]. During this process it is crucial to detect singularities along the way in order to recognize other merging branches of solutions which would otherwise be missed.

It is often useful to classify bifurcations into two principal classes: changes in the topology of the phase portrait [24] and changes of the solution set. For example, in the recent paper [25], the authors consider a reaction-diffusion model and show that by varying three model parameters simultaneously, they could delimit several bifurcation surfaces with high codimension such as transcritical, Bogdanov-Takens, Hopf and saddle-node bifurcations. Motivated by biological models, a classification of Bogdanov-Takens singularities is presented, and some properties of two singularities of high codimension is proven [21]. For specific examples of neural network systems, high codimensional bifurcations are studied in [22, 23].

This paper deals with bifurcations of solutions of differential equations and in particular with those whose solutions are critical points of a functional. These bifurcations can be classified according to catastrophe theory [26] and the first three bifurcation classes of this classification are the so-called $A$-series, $D$-series and $E$-series bifurcations. In fact, all occurring bifurcation phenomena can be classified and related to catastrophe theory [27, 28] for special instances, e.g. for systems of differential equations with Hamiltonian structure and boundary condition of Lagrangian type such as many physically relevant boundary conditions like periodic, Dirichlet or Neumann boundary conditions for second order ODEs. In comparison to ODEs, the detection of bifurcations of solutions of PDEs poses additional challenges including infinite dimensionality and (after discretisation) the need to solve large systems of equations. As a prototype of a PDE problem, we consider the semilinear Poisson equation

[TABLE]

for $u$ defined on $\Omega\subset\mathbb{R}^{d},\lambda\in\mathbb{R}^{k}$. Its parameter-dependent solutions can be regarded as stationary solutions of an associated reaction-diffusion equation with many applications in the physical sciences.

It seems that maps of singularities in the Banach space setting were first considered by Ambrosetti and Prodi [29] where they studied the situation of the fold map. Since then, fold and cusp maps have been studied extensively and many characterizations of fold and cusp maps have been given [30, 31, 32, 33, 34].

According to the Ambrosetti-Prodi theorem [29], the map $F(u)=\Delta u+f(u,\lambda)$ between appropriate functional spaces is a global fold map, provided certain hypotheses such as the convexity of the function $f$ with respect to its first component are satisfied. In [35], the authors show under mild conditions that convexity is indeed necessary. If the convexity condition of $f$ is not fulfilled, there exists a point with at least four preimages under $F$ and $F$ generically admits cusps among its critical points.

Following the fold and the cusp, the two subsequent singularities are the swallowtail and the butterfly whose maps were characterized by Ruf [36]. As an application, an elliptic boundary value problem with cubic nonlinearity, given by $\Delta u+\lambda u-u^{3}+h=0$ on a bounded, smooth domain $\Omega\subset\mathbb{R}^{d}$ with either Dirichlet or Neumann boundary conditions, is often considered. Here, the forcing term $h$ is a given Hölder continuous function with exponent $\alpha\in(0,1)$. This PDE, equipped with Neumann boundary conditions, has been studied in [37, 38, 36] in detail. In [37], Ruf showed that for certain parameter values $\lambda$ the elliptic equation has at most three solutions for arbitrary forcing terms $h$. In the subsequent paper [38], it is shown that for other parameter values $\lambda$ a secondary bifurcation occurs. In particular, at least five solutions can occur for certain forcing terms. Ruf [36] also gave a characterization of the solution geometry for the elliptic equation with cubic nonlinearity (and Neumann boundary conditions), while, independently, the same problem (with Dirichlet boundary conditions) was studied in [39]. Church and Timourian also derived abstract conditions for the global equivalence of a nonlinear mapping to the fold map [33] and the cusp map [34]. A nice historical overview about progress in singularity theory in differential equations is provided in [40].

The singularities mentioned above can be understood as singularities of functionals $G\colon E\to E$, where the equation $G(u)=0$ corresponds to the weak form of the PDE. In contrast, we investigate singularities of functionals $S\colon E\to\mathbb{R}$ where the weak formulation of the PDE is recovered as $\mathrm{D}S(u)=0$, where $\mathrm{D}$ denotes the Fréchet derivative. These singularities are related to classical catastrophe theory. We will mainly consider the important class of $A$-series singularities which can be defined informally as solutions at which the linearised equation has a one-dimensional kernel. However, due to the unified treatment of bifurcations in catastrophe theory, similar methods can be used to detect other high codimensional bifurcations in catastrophe theory such as $D$-series bifurcations. Figure 1 shows illustrations of the first three $A$-series singularities: the fold ($A_{2}$), the cusp ($A_{3}$) and the swallowtail ($A_{4}$) bifurcation. Further elementary singularities are discussed in [41]. Numerical methods to calculate bifurcation points of systems of equations are discussed in [42, 43, 44, 45, 46, 47, 18, 48]. Fink et al. [45] present a general mathematical framework for the numerical study of bifurcation phenomena associated with parameter-dependent equations. Singularities are directly computed as solutions of a minimally augmented defining system [46, 47]. Besides, singularities of functionals with an application to PDEs are discussed in [49]. A numerical experiment for detecting a high codimensional $A$-series bifurcation in a one-dimensional example can be found in [50, §4], for instance.

The two main approaches to bifurcation of solutions to PDEs which have been considered in the past are generalized Lyapunov-Schmidt reductions [43, 44, 18, 51] and topological methods in the calculus of variations which so far can only access relatively simple bifurcations. The Lyapunov-Schmidt reduction can be used to study solutions to nonlinear equations when the implicit function theorem cannot be applied, and allows the reduction of infinite-dimensional equations in Banach spaces to finite-dimensional equations. However, Lyapunov-Schmidt reductions do not use the variational structure. This motivates the development of a method for detecting high codimensional $A$-series bifurcations in high-dimensional or even infinite-dimensional phase spaces which makes use of the underlying variational structure. Based on the Infinite-Dimensional Splitting Lemma [52], we derive an augmented system for $A$-series singularities and apply it to PDEs by considering their variational formulations. For the numerical implementation, we use the most natural approach and discretise with a variational method. This procedure leads to a numerical method for the detection of high codimensional bifurcations in catastrophe theory which can be applied to a large class of parameter-dependent PDEs. In particular, our original contributions are:

The use of the Infinite-Dimensional Splitting Lemma. 2. 2.

A unified and simplified treatment of all $A$-series bifurcations. 3. 3.

The presentation in Banach spaces, i.e. our results apply both to variational PDEs and their discretisations. 4. 4.

Further simplification for the parameter-dependent semilinear Poisson equation (both in the continuous and discrete setting). 5. 5.

A unified treatment of the continuous problem and its discretisation.

The paper is organised as follows. In Section 2 we recall an observation by Golubitsky and Marsden [52] that catastrophe theory, i.e. the classification of the bifurcation behaviour of critical points of smooth, real-valued functions on (finite-dimensional) spaces, applies to smooth, nonlinear functionals on Banach spaces. We use this result to derive explicit bifurcation test equations (also known as determining equations or augmented systems) for all A-series singularities, expressed in terms of the derivatives of the original functional. Based on these derived augmented systems of equations, we propose a numerical method for finding high codimensional bifurcations in parameter-dependent PDEs. This numerical scheme is illustrated in Section 3 where we consider a Bratu-type boundary value problem as an of a second order PDE and detect its high codimensional bifurcations numerically. Finally, we conclude in Section 4.

2 Augmented systems for nonlinear functionals

2.1 The Splitting Lemma in Banach spaces

Let $E$ be a Banach space, let $U\subset E$ be an open neighbourhood of the origin and let $S\colon U\to\mathbb{R}$ be a smooth function with $S(0)=0$. We define the following two assumptions.

Assumption 2.1**.**

There exists an inner product $\langle\cdot,\cdot\rangle_{E}$ on $E$ and a Fredholm operator $T\colon E\to E$ of index 0 such that

[TABLE]

In 2.1 the symbol $\mathrm{D}$ denotes the Fréchet derivative. The second order derivative $D^{2}S(0)$ of $S$ at [math] is a symmetric bilinear form on $E\times E$. The index-0 Fredholm operator $T$ is symmetric such that $E=\ker T\oplus\mathrm{rg}\,T$, where $\ker T$ denotes the kernel of $T$ and $\mathrm{rg}\,T$ the range of $T$. We denote elements $u\in E$ by its components $u=(x,y)$ w.r.t. the splitting, i.e. $x\in\ker T$ and $y\in\mathrm{rg}\,T$.

Assumption 2.2**.**

There exists a partial gradient $\nabla_{y}S\colon U\to\mathrm{rg}\,T$ with $\nabla_{y}S(0)=0$ such that

[TABLE]

Theorem 2.1** (Infinite-Dimensional Splitting Lemma [52]).**

If 2.1 and 2.2 hold then there exists a fibred change of coordinates $(\bar{x},\bar{y})=(x,\eta(x,y))$ on $U$ fixing $(0,0)\in E$ with $\eta\colon U\to\mathrm{rg}\,T$ such that $\mathrm{D}_{y}\eta(0,0)=\mathrm{Id}$ and such that $S$ takes the form

[TABLE]

for a smooth function $r\colon\mathcal{O}\to\mathbb{R}$, where $\mathcal{O}$ is an open neighbourhood of $0\in\ker T$ and $r(0)=0$, $\mathrm{D}r(0)=0$ and $\mathrm{D}^{2}r(0)=0$.

In Theorem 2.1 $\mathrm{D}_{y}\eta$ denotes the Fréchet derivative of $\eta$ in the direction of $\mathrm{rg}\,T$. For a discussion of the setting and examples refer to [53]. The function $r$ is defined as follows: by the implicit function theorem, there exists a unique, smooth function $F\colon U_{\ker T}\to U_{\mathrm{rg}\,T}$ for open neighbourhoods $U_{\ker T}\subset\ker T$ and $U_{\mathrm{rg}\,T}\subset\mathrm{rg}\,T$ of the origin in $\ker T$ and $\mathrm{rg}\,T$, respectively, such that

[TABLE]

The map $r\colon U_{\ker T}\to\mathbb{R}$ is given by

[TABLE]

We see that critical points of $S$ correspond to critical points of $r$ which is defined on a finite-dimensional space. Singularity theory for $S$ thus reduces to ordinary, finite-dimensional catastrophe theory [54, 55]. More precisely, to determine the singularities of $S$ it suffices to determine the singularities of $r$.

2.2 Augmented systems for $A$-series singularities

$A$-series singularities for a given functional $S$ are defined as follows.

Definition 2.2** ($A$-series singularity).**

Let $E$ be a Banach space, let $U\subset E$ be a neighbourhood of the origin and let $S\colon U\to\mathbb{R}$ be a smooth function with $S(0)=0$. Assume that 2.1 and 2.2 hold and consider the function $r$ from (2). The function $S$ has a singularity of type $A_{n}$ at 0 if $\ker T$ (with $T$ as in 2.1) is one-dimensional and

[TABLE]

Remark 2.3**.**

The singularity $A_{2}$ is referred to as fold, $A_{3}$ as cusp, $A_{4}$ as swallowtail and $A_{5}$ as butterfly singularity.

Remark 2.4**.**

The Infinite-Dimensional Splitting Lemma (Theorem 2.1) allows us to borrow the notions of catastrophe theory and to define singularities of real-valued functionals $S\colon E\to\mathbb{R}$ fulfilling 2.1 and 2.2. In applications, $\mathrm{D}S(u)(v)=0\,\forall v\in E$ is the weak formulation of a PDE. In the literature the gradient structure is typically not exploited for this purpose. Instead the more general class of singularities of functionals $\tilde{S}\colon E\to\tilde{E}$ between two Banach spaces $E$ and $\tilde{E}$ is considered as in [29, 30, 31, 32, 33, 34, 36, 37, 38, 39]. The general problem class $\tilde{S}(u)=0$ contains the class of catastrophe problems $\mathrm{D}S(u)=0$. However, since the class of functions $E\to\tilde{E}$ is richer, this leads to a different notion of $A$-series singularities since stability properties do not coincide. To illustrate this point, let $\tilde{E}=E$. A map $\tilde{S}\colon E\to E$ with a singularity that is persistent under small perturbations $\tilde{S}+\epsilon\tilde{P}$ with functions $\tilde{P}\colon E\to E$ is not necessarily of the form $\tilde{S}=\mathrm{D}S$. On the other hand a map $S\colon E\to\mathbb{R}$ with a singularity that is persistent under small perturbations $S+\epsilon P$ with functions $P\colon E\to\mathbb{R}$ does not necessarily yield a map $\tilde{S}=\mathrm{D}S\colon E\to E$ with a singularity that is persistent under small perturbations with functions of the bigger problem class $\tilde{P}\colon E\to E$ since perturbations with $\tilde{P}\not=\mathrm{D}P$ are allowed. Let us illustrate the different notions of singularities on the cusp singularity. As proved by Whitney, the cusp map $f\colon\mathbb{R}^{2}\to\mathbb{R}^{2}$ with $f(z,w)=(z^{3}+zw,w)$ is stable. If we plot the $z$-component of solutions $(z,w)$ to the system $f(z,w)=(\mu_{1},\mu_{2})$ over the $(\mu_{1},\mu_{2})$-plane then we obtain the plot in the centre of Figure 1. On the other hand, consider the map $h\colon\mathbb{R}^{2}\to\mathbb{R}$ with $h(z,w)=\frac{1}{4}z^{4}+\frac{1}{2}w^{2}$. The map $h$ has a cusp singularity at $(0,0)$. Its universal unfolding in catastrophe theory is given by $h_{\mu}(z,w)=\frac{1}{4}z^{4}+\frac{1}{2}\mu_{2}z^{2}+\mu_{1}z+\frac{1}{2}w^{2}$. If we plot the $z$-component of solutions $(z,w)$ to the system $\nabla h_{\mu}(z,w)=(0,0)$ over the $(\mu_{1},\mu_{2})$-plane then we also obtain the plot in the centre of Figure 1. Despite the same visualisation, the cusp map $f$ is not to be confused with the map $\nabla h$. The map $f$ does not have a primitive $g\colon\mathbb{R}^{2}\to\mathbb{R}$ with $\nabla g=f$. Moreover, the map $f$ is stable in the class of smooth functions $\mathbb{R}^{2}\to\mathbb{R}^{2}$ while the cusp catastrophe $h\colon\mathbb{R}^{2}\to\mathbb{R}$ needs to be unfolded to be stable in the class of smooth functions $\mathbb{R}^{2}\to\mathbb{R}$. In this paper we will investigate the catastrophe setting and exploit the gradient structure for numerical continuation methods.

The $k^{\mathrm{th}}$ derivative of $S$ at a point $z\in U$ is a symmetric multi-linear form which we denote by $\mathrm{D}^{(k)}S(z)$. We can interpret the multi-linear form $\mathrm{D}^{(k)}S(z)$ as a linear form on $E^{\otimes k}=\bigotimes_{j=1}^{k}E$ where $\otimes$ denotes the tensor product. As a shorthand we define $S^{(k)}:=\mathrm{D}^{(k)}S$ and $S_{0}^{(k)}:=\mathrm{D}^{(k)}S(0)$.

We can express the condition for a function $S$ to have an $A$-series singularity in terms of (Fréchet)-derivatives of $S$. For this we define the multi-index set $\mathcal{J}_{k}^{n}$ for $n,k\in\mathbb{N}$ with $k\leq n$ as

[TABLE]

Moreover, for a multi-index $j\in\mathcal{J}_{k}^{n}$ we define $j!:=j_{1}!j_{2}!\cdots j_{n-k+1}!$.

Theorem 2.5**.**

Let $E$ be a Banach space and consider an inner product $\langle\cdot,\cdot\rangle_{E}$ on $E$. Let $U\subset E$ be a neighbourhood of the origin and let $S\colon U\to\mathbb{R}$ be a smooth function with $S(0)=0$. Consider the following algorithm consisting of a sequence of tests. The algorithm terminates and returns the current value of the integer $n$ if a test fails. In the algorithm $F$ is considered as a function symbol of an unknown smooth map $I\to E$, where $I$ is a small open interval in $\mathbb{R}$ containing 0 and with $F(0)=F^{(0)}(0)=0$, $F^{(1)}(0)=0$.

• •

Set $n=1$. Test $S_{0}^{(1)}=0$.

• •

Set $n=2$. Test whether the kernel of $a\mapsto S_{0}^{(2)}(a,\cdot)\in E^{\ast}$ is 1-dimensional.

• •

Set $n=3$. Select an element $\alpha\in E\setminus\{0\}$. Test $S^{(3)}(\alpha^{\otimes 3})=0$.

• •

Loop through the following two steps.

Set $n:=n+1$. Determine $F^{(n-2)}(0)\in\{{s}\alpha\,|\,{s}\in\mathbb{R}\}^{\bot_{E}}$ using

[TABLE]

[TABLE] 2. 2.

Test

[TABLE]

The following statements hold true.

• •

The algorithm returns $n=1$ if and only if [math] is not a critical point of $S$.

• •

The algorithm returns $n=2$ or does not terminate if and only if [math] is a critical point of $S$ but $S$ does not have an $A$-series singularity at [math].

• •

If 2.1 and 2.2 are satisfied w.r.t. $\langle\cdot,\cdot\rangle_{E}$ then the algorithm returns $n\geq 3$ if and only if $S$ has a singularity of type $A_{n-1}$.

Before proving the theorem, let us formulate some corollaries which illustrate how the conditions in Theorem 2.5 simplify for small values of $n$.

Definition 2.6**.**

We say that a functional $S$ fulfilling 2.1 and 2.2 has a singularity of type at least $A_{n}$ if

• •

it has a singularity of type $A_{N}$ with $N\geq n$ or

• •

the algorithm in Theorem 2.5 does not terminate.

Corollary 2.7** (Fold ($A_{2}$)).**

Let $E$ be a Banach space and $S$ be a (nonlinear) smooth functional defined on an open neighbourhood of $0\in E$. Assume that 2.1 and 2.2 hold for a Fredholm operator $T$ such that

[TABLE]

is satisfied for a non-trivial element $\alpha\in E$. Then $S$ has a singularity of type at least $A_{2}$ (fold singularity).

Corollary 2.8** (Cusp ($A_{3}$)).**

Let $E$ be a Banach space and $S$ be a smooth functional defined on an open neighbourhood of $0\in E$. Assume that 2.1 and 2.2 hold for a Fredholm operator $T$ and $\alpha\in E\setminus\{0\}$ such that (4) holds, i.e. $S$ has a singularity of type at least $A_{2}$. The functional $S$ has a singularity of type at least $A_{3}$ if and only if

[TABLE]

Corollary 2.9** (Swallowtail ($A_{4}$)).**

Let $E$ be a Banach space and $S$ be a smooth functional defined on an open neighbourhood of $0\in E$. Assume that 2.1 and 2.2 hold for a Fredholm operator $T$ and $\alpha\in E\setminus\{0\}$ such that (4) and (5) are satisfied, i.e. $S$ has a singularity of type at least $A_{3}$. The functional $S$ has a singularity of type at least $A_{4}$ if and only if

[TABLE]

is solvable for $v\in E$ and

[TABLE]

Corollary 2.10** (Butterfly ($A_{5}$)).**

Let $E$ be a Banach space and $S$ be a smooth functional defined on an open neighbourhood of $0\in E$. Assume that 2.1 and 2.2 hold for a Fredholm operator $T$ and $\alpha\in E\setminus\{0\}$ such that (4) and (5) are satisfied. Furthermore, assume that (6) holds for some $v\in E$ and (7) is satisfied, i.e. the functional $S$ has a singularity of type at least $A_{4}$. The functional $S$ has a singularity of type at least $A_{5}$ if and only if

[TABLE]

is solvable for $w\in E$ and

[TABLE]

Remark 2.11**.**

Notice that we do not need to require $v\in\mathrm{rg}\,T$ in (6) or $w\in\mathrm{rg}\,T$ in (8) since $S^{(2)}_{0}(w+t_{1}\alpha,v+t_{2}\alpha)=S_{0}^{(2)}(w,v)$ and $S^{(3)}_{0}(\alpha,v+t_{1}\alpha,\xi+t_{2}\alpha)=S^{(3)}_{0}(\alpha,v,\xi)$ for $\alpha\in E$ satisfying (4) and (5) such that the test equations are defined invariantly.

Remark 2.12**.**

The equations in the loop section of the algorithm presented in Theorem 2.5 can be obtained from the ${(n-2)^{\mathrm{th}}}$ and $n^{\mathrm{th}}$ complete exponential Bell polynomials as we will see from Lemma 2.15, Lemma 2.18 and Remark 2.19.

Definition 2.13** (Complete exponential Bell polynomial).**

The $n^{\mathrm{th}}$ complete exponential Bell polynomial is given as

[TABLE]

with the multi-index set $\mathcal{J}_{k}^{n}$ as in (3) and $j!=j_{1}!j_{2}!\ldots j_{n-k+1}!$ for $j\in\mathcal{J}_{k}^{n}$.

The first five complete exponential Bell polynomials are given by

[TABLE]

Remark 2.14**.**

Complete exponential Bell polynomials appear as coefficients in the following formal power series.

[TABLE]

Moreover, the $n^{\mathrm{th}}$ complete exponential Bell polynomial encodes information on the number of ways a set containing $n$ elements can be partitioned into non-empty, disjoint subsets. For example we can read off from

[TABLE]

that there is

• •

1 partition consisting of 4 sets of cardinality 1,

• •

6 partitions into 3 sets of which 2 have cardinality 1 and the other one has cardinality 2,

• •

4 partitions into 2 sets of which 1 has cardinality 1 and the other one has cardinality 3,

• •

3 partitions into 2 sets of cardinality 2,

• •

and 1 partition consisting of 1 set of cardinality 4.

(See [56, 57], for instance).

Let us prepare the proof of Theorem 2.5 with two lemmas.

Lemma 2.15**.**

Let $E$ be a Banach space, $U\subset E$ an open subset and $\alpha\in E$. Consider smooth functions $S\colon U\to\mathbb{R}$ and $F\colon I\to E$, where $I$ is an open interval $I\subset\mathbb{R}$ such that $r(s):=S(s\alpha+F(s))$ is defined on $I$. For $n\in\mathbb{N}$ we have

[TABLE]

On the right hand side of the equation multiplications are interpreted as tensor products. Moreover, the symbol “+” is replaced by “$+S^{(\mathrm{degree})}(s\alpha+F(s))$”, where $\mathrm{degree}$ is the count of factors in the tensor product to which $S^{(\mathrm{degree})}(s\alpha+F(s))$ is applied. In other words

[TABLE]

where $j=(j_{1},\ldots,j_{n-k+1})\in\mathcal{J}_{k}^{n}$ is defined in (3) and $j!=j_{1}!j_{2}!\ldots j_{n-k+1}!$ for a multi-index $j\in\mathcal{J}_{k}^{n}$.

Corollary 2.16**.**

In the setting of Lemma 2.15, if $0\in U$ and $S_{0}=0$, $S_{0}^{(1)}=0$, $S_{0}^{(2)}(\alpha,\xi)=0$ for all $\xi\in E$ and $F(0)=F^{\prime}(0)=0$ then the first five derivatives of $r$ evaluated at 0 are given by

[TABLE]

Proof of LABEL:{prop:DerivativesBellPoly}.

The statement is an extension of Faà di Bruno’s formula [58, §12.3] from $E=\mathbb{R}$ to arbitrary Banach spaces $E$. Since the chain rule for differentiation is valid for Fréchet derivatives [59, Thm. 2.47], the proof is analogous. ∎

As we see explicitly for $i\leq 5$ in Corollary 2.16, to determine $r^{(i)}(0)$ the $i-2$-jet of $F$ is required at 0 in the setting of Corollary 2.16. This holds in general.

Corollary 2.17**.**

In the setting of Corollary 2.16 the value $r^{(i)}(0)$ can be calculated from the $i-2$-jet of $F$ where $i\in\mathbb{N}$ with $i\geq 2$. We can write

[TABLE]

where $c_{1},c_{2}$ are any constants. On the right hand side of the equation multiplications are interpreted as tensor products. Moreover, the symbol “+” is replaced by “$+S^{(\mathrm{degree})}_{0}$”, where $\mathrm{degree}$ is the count of factors in the tensor product to which $S^{(\mathrm{degree})}_{0}$ is applied.

Proof.

We use the combinatorial interpretation of Bell polynomials and Lemma 2.15. If a partition of an $i$-set contains a subset with $i-1$ elements then there must be exactly one more non-empty subset containing one element to form a valid partition. Therefore, $F^{(i-1)}(0)=\left.\frac{\mathrm{d}^{i-1}}{\mathrm{d}s^{i-1}}\right|_{0}(s\alpha+F(s))$ only occurs together with $\alpha=\left.\frac{\mathrm{d}^{1}}{\mathrm{d}s^{1}}\right|_{0}(s\alpha+F(s))$ as $iF^{(i-1)}(0)\otimes\alpha$. This terms becomes an input argument of $S^{(2)}_{0}$ and, therefore, vanishes by the choice of $\alpha$. If a partition of an $i$-set contains a subset with $i$ elements then there cannot be another non-empty subset in the partition. Therefore, $F^{(i)}(0)=\left.\frac{\mathrm{d}^{i}}{\mathrm{d}s^{i}}\right|_{0}(s\alpha+F(s))$ becomes an input argument of $S^{(1)}_{0}$ which is zero. ∎

Lemma 2.18**.**

Let $E$ be a Banach space and let $U$ be an open neighbourhood of $0\in E$. Consider a smooth function $S\colon U\to\mathbb{R}$ with $S(0)=0$ such that 2.1 and 2.2 hold for a Fredholm operator $T$ with

[TABLE]

for a non-trivial element $\alpha\in E$. There exists an open interval $I$ containing 0 and a unique function $F\colon I\to\mathrm{rg}\,T$ with $F(0)=0$, $F^{\prime}(0)=0$ s.t.

[TABLE]

Moreover, all derivatives $F^{(n)}(0)$ with $n\geq 2$ can be obtained successively from

[TABLE]

*where $j=(j_{1},\ldots,j_{n-k+1})\in\mathcal{J}_{k}^{n}$ as defined in (3) and ${j!}=j_{1}!j_{2}!\ldots j_{n-k+1}!$. *

Remark 2.19**.**

Relation (12) can be expressed as

[TABLE]

where multiplications in the definition of the Bell polynomial (9) are interpreted as tensor products. Furthermore, the symbol “+” is replaced by “$+S^{(\mathrm{degree})}(s\alpha+F(s))\xi\otimes$” and parenthesis may be added around the input arguments of the form $S^{(\mathrm{degree})}(s\alpha+F(s))$.

Corollary 2.20**.**

The relations for $n=2,3,4,5$ of Lemma 2.18 read

[TABLE]

Proof of Lemma 2.18.

The operator $T$ defines an isomorphism on $\mathrm{rg}\,T$. Thus, the implicit function theorem applies to $\nabla_{y}S$ and together with (11) provides the existence and uniqueness of $F$ with $F(0)=0$ and $F^{\prime}(0)=0$ in analogy to the proof of Theorem 2.1 which can be found in [52]. Let $\xi\in\mathrm{rg}\,T$. Differentiating

[TABLE]

repeatedly with respect to $s$ we obtain the following relations.

[TABLE]

We encounter the same combinatorical relations as in Lemma 2.15 such that differentiating (13) $n$ times gives

[TABLE]

where multiplications in the equation above are interpreted as tensor products and “$+$” is replaced by “$+S^{(\mathrm{degree})}(s\alpha+F(s))\xi\otimes$. One may add parenthesis around the input arguments of the form $S^{(\mathrm{degree})}(s\alpha+F(s))$. In other words the relations are given by

[TABLE]

An evaluation at $s=0$ yields the claimed formula. In the $n^{\mathrm{th}}$ step of differentiation the term $F^{(n)}(0)$ only occurs as an input argument of $S^{2}_{0}$ and not elsewhere. Since the symmetric operator $T$ restricted to $\mathrm{rg}\,T$ is an isomorphism on $\mathrm{rg}\,T$ this successively determines all derivatives of $F$ at [math]. ∎

Proof of Theorem 2.5 and its corollaries.

The first two statements of the theorem follow by definition. Assume that 2.1 and 2.2 hold for an operator $T$ with 1-dimensional kernel. Let $\alpha\in\ker\,T\setminus\{0\}$. To analyse the singularities of $S$ it suffices to analyse the singularities of the function $r\colon U_{\ker}\to\mathbb{R}$ provided by Theorem 2.1. Using the identification $\mathbb{R}\xrightarrow{\sim}\ker T$, $s\mapsto s\alpha$, we identify $U_{\ker}$ with an open interval $I$ containing 0 and obtain $r\colon I\to\mathbb{R}$. The function $r$ has the form $r(s)=S(s\alpha+F(s))$ for a smooth function $F\colon I\to\mathbb{R}$. The functional $S$ has a singularity of type $A_{m}$ at 0 if and only if $r^{(k)}(0)=0$ for all $k\leq m$ and $r^{(m+1)}(0)\not=0$.

The algorithm presented in the statement of Theorem 2.5 consists of a sequence of tests and a variable $n$ acts as a counter. If the state of the variable $n$ is $k$ then the test in the algorithm corresponds to testing $r^{(k)}(0)=0$. This can be seen from Lemma 2.15. (The formula for $r^{(k)}(0)$ is related to the $k^{\mathrm{th}}$ complete exponential Bell polynomial.) To evaluate $r^{(k)}(0)$ the $k-2$-jet of $F$ is required as observed in Corollary 2.17. If the state of the variable $n$ is $k\geq 4$ then $F^{(k-2)}(0)$ gets determined in the algorithm just before $r^{(k)}(0)=0$ is tested. Lemma 2.18 justifies that the algorithm can determine $F^{(k-2)}(0)$ via the given formula (which is related to the $(k-2)^{\mathrm{th}}$ complete exponential Bell polynomial) if all values $F^{(i)}(0)$ are defined for $i\leq k-3$. The values $F(0)$ and $F^{\prime}(0)$ are set to be 0. The theorem follows by induction.

Corollaries 2.7, 2.8, 2.9 and 2.10 follow from Corollaries 2.16 and 2.20. In Corollary 2.9 and Corollary 2.10 the bifurcation test equations $r^{(4)}(0)=0$ and $r^{(5)}(0)=0$ have been simplified using (6) and (8). ∎

Proposition 2.21**.**

Let $E$ be a Banach space and $S$ be a real-valued functional defined on an open neighbourhood of $0\in E$ such that 2.1 and 2.2 hold. Consider the function $r$ of the proof of Theorem 2.5, whose derivatives are bifurcation test equations. If $S$ has a singularity of type $A_{2k+1}$ with $k\in\mathbb{N}$ at [math] then the signature of $r^{(2n+2)}(0)$ is well-defined.

Proof.

The statement follows from the classification of singularities up to right-equivalence in catastrophe theory [60] or can be deduced from our considerations as follows. The map $r(s)=S(s\alpha+F(s))$ from the proof of Theorem 2.5 is defined uniquely up to the choice of $\alpha$, where $\alpha$ is as in 2.1. The determining equations for the jet of $F$ at 0 and the bifurcation test equations $r^{(j)}(0)$ ($j\in\mathbb{N}$) are related to Bell polynomials by Lemmas 2.15, 2.18 and 2.19.

• •

In any partition of an even amount of elements there must be an even number (or none) of subsets with odd cardinality.

• •

In any partition of an odd amount of elements there must be an odd number of subsets with odd cardinality.

Using the two combinatorial observations above we see inductively that the signatures of the derivatives $F^{(2j)}(0)$ are defined invariantly of $\alpha$, all derivatives $F^{(2j+1)}(0)$ change to $-F^{(2j+1)}(0)$ ($j\in\mathbb{N}$) as $\alpha\mapsto-\alpha$ and can conclude that the signature of $r^{(2k+2)}(0)$ is well-defined. ∎

Remark 2.22**.**

The signature considered in 2.21 occurs in catastrophe theory if a classification of singularities up to right-equivalence is considered [60]. The singularities $A_{2k+1}$ do not have a signature. If the algorithm in Theorem 2.5 returns $n=2k+2$ with $k\in\mathbb{N}$ and 2.1 and 2.2 are satisfied then the singularity $A_{2k+1}$ is of the positive type if and only if the last test equation (which corresponds to the $n^{\mathrm{th}}$ complete exponential Bell polynomial) is positive. Otherwise the singularity is of the negative type.

Remark 2.23**.**

If the functional $S$ has $2k+1$ parameters ($k\in\mathbb{N}$) then, under non-degeneracy conditions, $A_{2k+1}$ singularities occur as 1-parameter families (by the implicit function theorem). If two branches of $A_{2k+1}$ singularities merge in a $A_{2k+2}$ singularity then one consists of singularities of the positive type and the other one of singularities of the negative type. This is because the bifurcation test equation $r^{(2k+2)}$, which determines the signs of the $A_{2k+1}$ singularities, must have a non-degenerate zero at the $A_{2k+2}$ singularity, i.e. its graph intersects the axis of abscissas transversally.

Remark 2.24**.**

Remark 2.23 applies to the fold bifurcation $A_{2}$ as well. In the finite-dimensional case the signature of a solution $z$ can be obtained as the sign of the determinant222which does not depend on the choice of basis since $\mathrm{D}S(z)=0$ of the Hessian matrix $S^{(2)}_{z}$ of $S$ at $z$. In a numerical computation the sign can be determined by performing an LU-decomposition of $S^{(2)}_{z}$ without pivoting and counting whether the number of positive signs on the diagonal of $U$ is even or odd. Keeping track of the signatures of solutions, cusps,…, provides information on which ones may be able to meet in a bifurcation.

3 Example: semilinear Poisson equation

We will exemplify how the augmented systems derived in Section 2.2 can be applied to PDEs. For this, we consider a second order, semilinear PDE describing the steady state solutions in a reaction-diffusion process [61]. First, we will justify that the theory presented in Section 2.1 applies and write down the continuous recognition equations. By a concrete numerical example we will show how augmented systems can be employed in continuation methods to find high codimensional singularities.

3.1 Setup

For a smooth function $f\colon\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}$ we consider the homogeneous Dirichlet problem

[TABLE]

on an open and bounded domain $\Omega\subset\mathbb{R}^{d}$ with boundary $\partial\Omega$ of class $\mathcal{C}^{k}$, where $k>\frac{d}{2}+1$. We denote the standard volume form on $\Omega$ by $\mathrm{d}{\bf x}$. In analogy to [53, Example 7] we consider the following setting which will allow us to employ the Infinite-Dimensional Splitting Lemma. See [36], for instance, for an alternative treatment. The Sobolev space $H^{k}(\Omega)$ is compactly embedded into $\mathcal{C}^{1}(\overline{\Omega})$ [62, Thm 6.2]. Consider the Hilbert space $E=H_{0}^{1}(\Omega)\cap H^{k}(\Omega)$ with the structure inherited from $H^{k}(\Omega)$ [62, 63]. Let $\bar{f}$ be s.t. $\frac{\partial}{\partial t}\bar{f}(t,\lambda)=f(t,\lambda)$ and consider the non-linear functional $S\colon E\times\mathbb{R}^{s}\to\mathbb{R}$ defined as

[TABLE]

The Fréchet derivatives of $S$ in the directions $v_{1},v_{2},\ldots\in E$ exist and are given as

[TABLE]

Here $f^{(l)}(t,\lambda)=\frac{\partial^{l}}{\partial t^{l}}f(t,\lambda)$. The equation $S^{(1)}(u,\lambda)(v)=0$ for all $v\in E$ is a weak formulation of (14).

We consider the bilinear form $\langle\cdot,\cdot\rangle_{E}\colon E\times E\to\mathbb{R}$ with

[TABLE]

where $\nabla v$ and $\nabla w$ denote weak derivatives of $v,w\in E$ and $\langle\cdot,\cdot\rangle$ the scalar product in $\mathbb{R}^{d}$. The bilinear form $\langle\cdot,\cdot\rangle_{E}$ is symmetric, positive definite by Poincaré’s inequality and bounded using the Cauchy-Schwarz inequality. The embedding $E=H_{0}^{1}(\Omega)\cap H^{k}(\Omega)\hookrightarrow H^{k-2}(\Omega)$ is compact [62, Thm 6.2]. Moreover, the Dirichlet Laplacian $L$ is an isomorphism $H_{0}^{1}(\Omega)\cap H^{k}(\Omega)\to H^{k-2}(\Omega)$ [64]. Therefore, the operator $\Delta^{-1}\colon E\to E$ defined as the composition

[TABLE]

is compact.

For each $(u,\lambda)\in E\times\mathbb{R}^{s}$ define the operator $T(u,\lambda)\colon E\to E$ by

[TABLE]

For each $(u,\lambda)\in E\times\mathbb{R}^{s}$ the operator $T{(u,\lambda)\colon E\to E}$ is a Fredholm operator: since $v\mapsto f^{\prime}(u,\lambda)v$ is continuous from $E$ into $E$ and $\Delta^{-1}\colon E\to E$ is compact, it follows that the composition is compact and $T(u,\lambda)\colon E\to E$ is a Fredholm operator of index 0 (Fredholm alternative). We have

[TABLE]

The operator $T(u,\lambda)$ is symmetric w.r.t. $\langle\cdot,\cdot\rangle_{E}$ such that we obtain the $\langle\cdot,\cdot\rangle_{E}$-orthogonality of $\ker(T(u,\lambda))$ and $\mathrm{rg}\,(T(u,\lambda))$. Since $T(u,\lambda)$ is a Fredholm operator, both spaces are closed in $E$. Moreover, $\mathrm{rg}\,(T(u,\lambda))$ has finite codimension. Since $T(u,\lambda)$ is symmetric, it is an elementary exercise to deduce that $E=\ker(T(u,\lambda))\oplus\mathrm{rg}\,(T(u,\lambda))$. The projection $\mathrm{pr}\colon E\to\mathrm{rg}\,(T(u,\lambda))$ induced by the splitting is continuous. For each $\lambda\in\mathbb{R}^{s}$ we define the operator $\nabla_{y}S(\lambda)\colon E\to\mathrm{rg}\,T(u,\lambda)$ as

[TABLE]

to $\mathrm{rg}\,T(u,\lambda)$. We have

[TABLE]

for all $v\in\mathrm{rg}\,(T(u,\lambda))$. The observations imply the following proposition.

Proposition 3.1**.**

The equation

[TABLE]

is a weak formulation of (14). Furthermore, assuming that $S^{(1)}(u,\lambda)=0$, 2.1 and 2.2 hold with $\langle\cdot,\cdot\rangle_{E}$ (15), the operator $T{(u,\lambda)\colon E\to E}$ and $\nabla_{y}S{(\lambda)\colon E\to\mathrm{rg}\,T(u,\lambda)}$ defined above for each $(u,\lambda)\in E\times\mathbb{R}^{s}$.

We conclude that Theorem 2.1 applies to $S$ such that the bifurcation behaviour of (14) reduces to finite-dimensional catastrophe theory.

3.2 Augmented systems for the example problem

Let us write down the augmented systems for (14) provided by Corollaries 2.7, 2.8, 2.9 and 2.10.

Solution

[TABLE]

Fold $(A_{2})$

[TABLE]

Cusp $(A_{3})$

[TABLE]

Swallowtail $(A_{4})$

[TABLE]

[TABLE]

Butterfly $(A_{5})$

[TABLE]

[TABLE]

Remark 3.2**.**

We can impose

[TABLE]

in (19) or $\langle v,w\rangle_{E}=0$ in (20) as discussed in Remark 2.11. The uniqueness condition allows us to interpret $v$ as $F^{\prime\prime}(0)$ and $w$ as $F^{\prime\prime\prime}(0)$.

Singularity $A_{n}$, $n\geq 4$

The equation to determine $F^{(n-2)}(0)\colon\Omega\to\mathbb{R}$ reads

[TABLE]

To obtain the correct expression for the term $B_{n-2}(\alpha,F^{\prime\prime}(0),\ldots,F^{(n-2)}(0))$ in (21), the Bell polynomial $B_{n-2}$ is first presented in its monomial form as in (10), then the summation sign $+$ is replaced by $+f^{(\mathrm{degree})}(u,\lambda)$, where $\mathrm{degree}$ denotes the degree of the monomial it is multiplied with. Finally, the arguments $\alpha,F^{\prime\prime}(0),\ldots,F^{(n-2)}(0)$ are substituted into the expression. The bifurcation test equation is given as

[TABLE]

(cf. Corollary 2.17)

To obtain the correct expression for the term $B_{n}(\alpha,F^{\prime\prime}(0),\ldots,F^{(n-2)}(0),0,0)$ in (22) the Bell polynomial $B_{n}$ is first presented in its monomial form. Summands, which consist of a degree 2 monomial $cx_{i}x_{j}$ are replaced by

[TABLE]

where $c$ denotes the occurring factor. In the other summands we add $f^{(\mathrm{degree}-1)}(u,\lambda)$ as a factor, whereas $\mathrm{degree}$ is the degree of the monomial making up the summand.

Remark 3.3**.**

We see that the fact that $S$ depends only quadratically on $u$ has led to a significant simplification compared to the formulas for general functionals. Indeed, the first equation in (21) is of the form

[TABLE]

for a smooth map $g\colon\Omega\to\mathbb{R}$ and with $\kappa(x)=f^{\prime}(u,\lambda)$. The map $F^{(n-2)}(0)\colon\Omega\to\mathbb{R}$ is sought. While the initial PDE (14) is a semilinear Poisson equation [65, 66], the equations which are added in the augmented system are linear PDEs. If $\kappa(x)\geq 0$ or $\kappa(x)\leq 0$ for all $x\in\Omega$ then this is the generalised Poisson equation considered in [67].

Example 3.4**.**

Let $f(u,\lambda)=\sum_{l=1}^{n}\frac{1}{l!}\lambda_{l}u^{l}$. For any choice of $\lambda$ the constant function $u=0$ is a solution to the problem (14). The condition for a singularity of type $A_{n}$ at $(\lambda,u)=(\lambda,0)$ is fulfilled if and only if $-\lambda_{1}$ is a simple eigenvalue of the Dirichlet Laplacian $\Delta\colon H_{0}^{1}(\Omega)\cap H^{k}(\Omega)\to H^{k-2}(\Omega)$ and

[TABLE]

provided that $\int_{\Omega}\alpha^{i}\mathrm{d}{\bf x}\not=0$ for $3\leq i\leq n+1$, where $\alpha$ is the eigenfunction to the eigenvalue $-\lambda_{1}$.

Proof.

The function $u=0$ solves (16). We have $f^{(l)}(0,\lambda)=\lambda_{l}$. The system (17) reads

[TABLE]

It has a unique solution if and only if $-\lambda_{1}$ is a simple eigenvalue of the Dirichlet Laplacian. The cusp condition simplifies to $\lambda_{2}=0$. Assuming (16), (17), (18) we see that $F^{\prime\prime}(0)=0$ solves (21) with $n=4$. The swallowtail condition is fulfilled if and only if $\lambda_{3}=0$. Assuming that $F^{\prime\prime}(0)=\ldots=F^{(t-2)}(0)=0$ and $\lambda_{3}=\lambda_{4}=\ldots=\lambda_{t-1}=0$ we see that $F^{(t-1)}(0)=0$ solves (21) with $n=t+1$. The bifurcation test equation (22) for $A_{t+1}$ is fulfilled if and only if $\lambda_{t}=0$. The claim follows by induction. ∎

3.3 Numerical experiment

The classical Bratu problem considers the PDE

[TABLE]

on a $d$-dimensional cube with zero Dirichlet boundary values. The boundary value problem is popular to study fold and cusp bifurcations [68]. In the following we consider the domain $\Omega=(0,1)\times(0,1)$ and set

[TABLE]

in the Dirichlet problem (14) given as

[TABLE]

The considered problem coincides with the Bratu problem for $\lambda_{3}=0$. Numerical experiments have been performed in [69]. The third parameter has been added to create a swallowtail bifurcation, which we will find numerically. As the topological boundary $\partial\Omega$ of the domain is not regular, the assumptions in Section 3.1 do not hold. However, the following numerical experiment illustrates on a classical example how the derived augmented systems can be used to locate bifurcations. As before, the primitive w.r.t. the first argument is denoted by $\bar{f}$, i.e. $\frac{\partial}{\partial t}\bar{f}(t,\lambda)=f(t,\lambda)$.

3.3.1 Discretisation via a discrete Lagrangian method

In view of Sections 2 and 3.1 it appears natural to use a variational method (see [70, VI.6] or [71], for instance) to discretise (14). We obtain a mesh on $\Omega$ as the cross product of a uniform mesh in the $x$ direction with $N$ interior mesh points and spacing $\Delta x=\frac{1}{N+1}$ and a uniform mesh in the $y$ direction with $M$ interior mesh points and spacing $\Delta y=\frac{1}{M+1}$. Real-valued, continuous functions $u\colon\Omega\to\mathbb{R}$ are represented as matrices $U\in\mathbb{R}^{N\times M}$ where the component $U_{i,j}$ corresponds to the value $u(i\Delta x,j\Delta y)$ on the interior grid. Alternatively, we can flatten the matrix $U$ to a vector $\bar{u}\in\mathbb{R}^{N\cdot M}$ such that $U_{i,j}=\bar{u}_{(j-1)N+i}$. For the functional

[TABLE]

from Section 3.1 we consider the discrete functional

[TABLE]

with $U_{0,j}=U_{i,0}=U_{N+1,j}=U_{i,M+1}=0$. For $k\in\mathbb{N}$ define

[TABLE]

The matrices

[TABLE]

are the standard central finite-difference discretisations of the operators $\frac{\partial^{2}}{\partial x^{2}}$ and $\frac{\partial^{2}}{\partial y^{2}}$. Define

[TABLE]

where $\otimes$ denotes the Kronecker tensor product. Notice that

[TABLE]

In the expression above the function $\bar{f}$ is evaluated component-wise, i.e the $k^{\text{th}}$ component of $\bar{f}(u,\lambda)$ is given by $\bar{f}(u_{k},\lambda)$. We have

[TABLE]

Remark 3.5**.**

The expression (23) is the flattening of the matrix-valued function

[TABLE]

which can be obtained by a direct discretisation of the PDE (14) using standard central difference approximations for second derivatives. We see that the equation $\mathrm{D}_{u}S_{\Delta}(\bar{u},\lambda)=0$ obtained from an application of the discrete Lagrangian principal to $S$ leads to the same equations as a discretisation of the Laplacian operator using central finite-differences.

3.3.2 Application of pseudoarclength continuation to augmented discrete systems

To locate a singularity of a given type one can write down the augmented system of the singularity and solve the system using an iterative solver. This is called a direct method. For the convergence of the solver a good initial guess is required which is usually not available a priori. Instead, one may employ a continuation method and first continue a known solution along one parameter until a fold singularity is detected, then augment the system using the fold bifurcation test equation and continue a line of folds until a cusp is detected and so on. In the following numerical experiment we will use pseudoarclength continuation [72]. Our strategy of performing conceptually simple one-dimensional continuations of singularities of high codimension can be contrasted to higher-dimensional continuation of solutions [73]. More information on continuation methods can be found in [74, 69, 75, 5].

Whenever a bifurcation is detected one has the option to locate the singularity of the discrete system exactly using a direct method. Generally speaking, unless one has prior knowledge about the bifurcation diagram or is interested only in a specific parameter region one needs to search in different directions for bifurcations. Once arrived at a high codimensional bifurcation, the discretisation parameter is decreased and a direct method is applied to approximate the location of the singularity in the continuous system. As starting values interpolated data from the coarser systems can be utilised.

To simplify notation in the following we will write $u$ instead of $\bar{u}$, $G$ for $\mathrm{D}S_{\Delta}$ and $G_{u}$ for the Jacobian matrix $\mathrm{D}_{u}G$.

Fold

Fixing $(\lambda_{2},\lambda_{3})=(0,0)$ and starting at $\lambda_{1}=0$ we continue the solution $u=0\in\mathbb{R}^{N\cdot M}$ by applying pseudoarclength continuation to

[TABLE]

until we find a fold. The singularity is detected when $\lambda_{1}(s)$ changes from being increasing to decreasing in $s$. (See the left plot in Figure 2.)

Cusp

We allow $\lambda_{2}$ to be parameter dependent as well and apply pseudoarclength continuation to

[TABLE]

starting with the approximated fold data and a random, normalised guess for $\alpha$. (See the right plot in Figure 2.) During the continuation process we monitor the cusp condition

[TABLE]

where raising $\alpha$ to the third power is to be understood component-wise. After detecting a change of sign in the left-hand side of (24), we improve the accuracy of the location of the singular point by solving

[TABLE]

using Newton iterations.

Swallowtail

We allow an $s$-dependence of $\lambda_{3}$ and apply pseudoarclength continuation to the system

[TABLE]

starting from the calculated cusp position. During this process, we monitor the swallowtail condition

[TABLE]

where $f_{uu}(u,\lambda).\alpha^{2}$ denotes the component-wise product of the vectors $f_{uu}(u,\lambda)$ and $\alpha^{2}$. In each continuation step the vector $v$ is obtained as follows: we calculate the vector $(v^{T},t)^{T}\in\mathbb{R}^{N\cdot M+1}$ which minimises the euclidean norm of

[TABLE]

For this we calculate the (under non-degeneracy conditions unique) solution to

[TABLE]

and obtain $v$ as $v=G_{u}(u,\lambda_{1},\lambda_{2},\lambda_{3})^{T}\bar{v}$. (Notice that transposition can be neglected since $G_{u}(u,\lambda_{1},\lambda_{2},\lambda_{3})$ is symmetric.)

As the left-hand side of (25) changes sign, we detect a candidate for a swallowtail point at $(u_{\mathrm{SW}},\lambda^{\mathrm{SW}})=(u_{\mathrm{SW}},\lambda^{\mathrm{SW}}_{1},\lambda^{\mathrm{SW}}_{2},\lambda^{\mathrm{SW}}_{3})$. (See Figure 5.) Indeed, the swallowtail condition (25) has a regular root at the swallowtail point which implies that the singularity is not further degenerate.

In Figure 6 we do a fold continuation using pseudoarclength continuation applied to $F^{(1)}$ starting at a cusp point near $(u_{\mathrm{SW}},\lambda^{\mathrm{SW}})$ while fixing $\lambda_{3}$. The fold line is continued in both directions and shows the characteristics of a swallowtail bifurcation. This verifies that $(u_{\mathrm{SW}},\lambda^{\mathrm{SW}})$ is indeed a swallowtail point of the discretised system.

Remark 3.6**.**

As is known from catastrophe theory [76, 60], the only generic333persistent under small perturbations. For exact notions see [76], for instance. bifurcations of codimension smaller or equal to 3 which critical points of a function $R\colon\mathbb{R}^{n}\to\mathbb{R}$ can undergo are $A_{2}$, $A_{3}$, $A_{4}$, $D^{+}_{4}$ and $D^{-}_{4}$. However, if a singularity $D^{+}_{4}$ or $D^{-}_{4}$ occurs along a line of cusp bifurcations then the swallowtail condition tends to $+\infty$ or $-\infty$ as one approaches the singularity. (See Figures 3 and 4.) Moreover, at the singularity, $\mathrm{Hess}(R)(z)$ has a 2-dimensional kernel and the value of the swallowtail condition is not defined.

Remark 3.7**.**

Forming augmented systems for the discrete functional $S_{\Delta}$ (recall $\mathrm{D}_{u}S_{\Delta}=G$) gives rise to the same system of equations as discretising the continuous augmented systems from Section 3.2. In other words, forming the augmented system commutes with discretisation. This extends Remark 3.5.

3.3.3 Determination of the position and convergence of the swallowtail point

To calculate the position of the bifurcation point more accurately we apply Newton iterations to $F^{(3)}(u,\alpha,\bar{v},\lambda)$ given as

[TABLE]

until convergence using $(u_{\mathrm{SW}},\alpha_{\mathrm{SW}},0_{NM\times 1},\lambda^{\mathrm{SW}})$ as initial guess. We obtain a root $(u^{1}_{\mathrm{SW}},\alpha^{1}_{\mathrm{SW}},v^{1}_{\mathrm{SW}},\lambda^{1}_{\mathrm{SW}})$ with $(u^{1}_{\mathrm{SW}},\alpha^{1}_{\mathrm{SW}},\lambda^{1}_{\mathrm{SW}})$ near $(u_{\mathrm{SW}},\alpha_{\mathrm{SW}},\lambda^{\mathrm{SW}})$. We successively increase $(N,M)=(N,N)$ and repeat the process of finding a root of $F^{(3)}$. As initial guesses for Newton’s method in the $j^{\mathrm{th}}$-step we use $\lambda^{j}_{\mathrm{SW}}$ from the previous calculation and linear interpolations of $u^{j}_{\mathrm{SW}}$, $\alpha^{j}_{\mathrm{SW}}$ and $v^{j}_{\mathrm{SW}}$ on the new grid with Dirichlet boundary conditions. The Jacobian matrix of $F^{(3)}$ required for the Newton iterations is given as

[TABLE]

with

[TABLE]

Here we use the convention that “.” denotes point-wise multiplication. Point-wise multiplication of a column vector with a matrix means that the vector is multiplied point-wise with each column of the matrix. Analogously for row vectors. Moreover, the power of a vector is to be understood point-wise. Zero matrices of dimension $s\times t$ are denoted by $0_{s\times t}$ and identity matrices by $\mathrm{Id}_{s\times t}$.

Remark 3.8**.**

Appropriate sub-matrices of $\mathrm{D}F^{(3)}$ correspond to $\mathrm{D}F^{(1)}$ and $\mathrm{D}F^{(2)}$. These have been used for the pseudoarclength continuation described in Section 3.3.2. The Jacobian matrices $\mathrm{D}F^{(1)}$, $\mathrm{D}F^{(2)}$ and $\mathrm{D}F^{(3)}$ and the matrices $L$, $G_{u}$ involved in the function evaluations $F^{(1)}$, $F^{(2)}$ and $F^{(3)}$ are sparse and represented using an appropriate datatype in the numerical calculations. Moreover, compared to equations for a general functional $S$, the structure of the semilinear Poisson equation dramatically simplifies occurring equations and numerical complexity because the multilinear operators $G_{uu}$, $G_{uuu}$, …are diagonal and their diagonal is given by an evaluation of an appropriate derivative of $f$ at $u$. The simplicity of the structure has been observed in the continuous setting in Remark 3.3.

Figure 7 shows that as $N=M$ increases, the position of the swallowtail point converges in the parameter space. We also observe that the data $(u,\alpha,v)$ converges to a fixed shape. Figure 9 shows the shape of $(u,\alpha,v)$ at the approximated swallowtail point for $(N,M)=(85,85)$.

4 Conclusion

In conclusion, we derived bifurcation test equations for $A$-series singularities of nonlinear functionals and, based on these equations, we developed numerical methods for the detection of high codimensional branching bifurcations in parameter-dependent PDEs where a variational integrator is used for the discretisation of the problem. This numerical computation is illustrated by detecting a swallowtail bifurcation in a Bratu-type problem. As part of future research, numerical experiments for finding high codimensional bifurcations for other PDE types can be performed. Another interesting question is the rigorous proof of the convergence rate to the singularity of the numerical method. Besides, bifurcation test equations for $D$-series singularities can be derived and the bifurcations can be numerically detected for different PDE types by considering either variational or non-variational integrators.

We would like to thank Wolf-Jürgen Beyn and Tom ter Elst for useful feedback on the manuscript and Hinke Osinga, Bernd Krauskopf, and Chris Budd for fruitful discussions. This research was supported by the Marsden Fund of the Royal Society Te Apārangi. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Geometry, compatibility and structure preservation in computational differential equations (2019) where work on this paper was undertaken. This work was supported by EPSRC grant No. EP/K032208/1. LMK was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L016516/1, the RISE projects CHiPS and NoMADS, and the German National Academic Foundation (Studienstiftung des Deutschen Volkes).

ORCID iDs

L M Kreusser https://orcid.org/0000-0002-1131-1125

R I McLachlan https://orcid.org/0000-0003-0392-4957

C Offen https://orcid.org/0000-0002-5940-8057