Numerical Continuation and SPDE Stabilityfor the 2D Cubic-Quintic Allen-Cahn Equation

Numerical Continuation and SPDE Stability
for the 2D Cubic-Quintic Allen-Cahn Equation

Christian Kuehn Institute for Analysis and Scientific Computing, Vienna University of Technology, 1040 Vienna, Austria

We study the Allen-Cahn equation with a cubic-quintic nonlinear term and a -trace-class stochastic forcing in two spatial dimensions. This stochastic partial differential equation (SPDE) is used as a test case to understand, how numerical continuation methods can be carried over to the SPDE setting. First, we compute the deterministic bifurcation diagram for the PDE, i.e. without stochastic forcing. In this case, two locally asymptotically stable steady state solution branches exist upon variation of the linear damping term. Then we consider the Lyapunov operator equation for the locally linearized system around steady states for the SPDE. We discretize the full SPDE using a combination of finite-differences and spectral noise approximation obtaining a finite-dimensional system of stochastic ordinary differential equations (SODEs). The large system of SODEs is used to approximate the Lyapunov operator equation via covariance matrices. The covariance matrices are numerically continued along the two bifurcation branches. We show that we can quantify the stochastic fluctuations along the branches. We also demonstrate scaling laws near branch and fold bifurcation points. Furthermore, we perform computational tests to show that, even with a sub-optimal computational setup, we can quantify the subexponential-timescale fluctuations near the deterministic steady states upon stochastic forcing on a standard desktop computer setup. Hence, the proposed method for numerical continuation of SPDEs has the potential to allow for rapid parametric uncertainty quantification of spatio-temporal stochastic systems.

Keywords: Allen-Cahn, SPDE, cubic-quintic Ginzburg-Landau, Lyapunov equation, numerical continuation, spectral noise approximation, uncertainty quantification, finite differences, iterative linear solvers, critical transitions, bifurcation diagram.

1 Introduction

Spatio-temporal systems with stochastic forcing and parameter dependence arise extremely frequently in practical applications. Stochastic partial differential equations (SPDEs) have emerged as one of the main examples for spatio-temporal stochastic systems, e.g. in ecological invasion waves [25, 50, 58], voter models [59, 20], surface and thin film growth models [17, 14], statistical mechanics [43, 37], fluid dynamics [41, 10], neuroscience [73, 67], and many general pattern formation problems [13, 29]. In this paper, we study as a prototypical model a version of the Allen-Cahn SPDE, which can formally be written as


where , , are parameters, are suitable maps, is the Laplacian and is a stochastic process. Our main interest is to use (1) as a test problem to illustrate a different approach to quantify stochastic fluctuations using numerical continuation methods. A more precise formulation of (1) and the detailed reasons for choosing a certain Allen-Cahn-type SPDE are discussed below.

Currently, a standard approach to study the dynamics of (1) numerically is centered around Monte-Carlo (MC) algorithms and their variations. The strategy is to use forward integration of the SPDE. This requires three different components:

  1. Sample over a sufficiently large number of initial conditions and determine which initial conditions lead to similar dynamics.

  2. Sample over a sufficiently large number of noise realizations for the process and compute the required averaged quantities of interest over several sample paths.

  3. Sample over the parameters and plot the bifurcation diagram at each point by repeating steps (MC1)-(MC2).

There are several strategies that partially alleviate some of the computational time required for carrying out (MC1)-(MC3), such as quasi-Monte-Carlo methods [18, 61] based on discrepancy theory or multilevel-Monte-Carlo [39, 31] based upon a coarse-/fine-sampling hierarchy. However, there are some drawbacks when following an MC-based strategy. Due to the three sampling spaces (initial conditions, noise, parameters), the dimension of the computational problem can become quite restrictive, particularly for the case of SPDEs. Furthermore, inspection of individual sample paths is often necessary to study dynamical objects such as bifurcations, metstable paths, or locally invariant manifolds. This can make MC methods in practical applications surprisingly labor-intensive. Even for deterministic differential equations, such as ordinary differential equations (ODEs) or partial differential equations (PDEs), this issue has been known for some time as discussed in [72, p.1-2].

Numerical continuation for ODEs and PDEs provides a systematic alternative to study dynamical phenomena. The idea is to track dynamical objects based upon parameter variation by a predictor-corrector scheme [44, 3, 34]. This approach is used successfully in a wide variety of applications to study the parameter depedence of models [47, 52, 70]. In combination with singularity detection algorithms, a common application of numerical continuation methods is to compute bifurcation diagrams [23, 53]. For ODE bifurcation analysis, excellent software packages are available such as AUTO [24] or MatCont [22]. Although there are many other applications of numerical continuation methods (see e.g. the short summary in the introduction to [51]), there is relatively little work, which tries to apply these ideas to stochastic problems [6, 48].

A possible strategy to make numerical continuation methods useful for stochastic ordinary differential equations (SODEs) suggested in [48] is composed of the following main steps:

  1. Continue the deterministic problem without forcing with techniques for ODEs. For example, one might want to compute the steady states and their bifurcations depending upon parameters.

  2. Use stochastic analysis on the level of the SODE to derive deterministic algebraic and/or differential equations for important stochastic quantities. For example, one may derive an ODE for the evolution of local stochastic fluctuations around a steady state for the linearized system. Another example are algebraic equations for the pre-factor and exponent in Kramers’ law governing metastable gradient systems.

  3. Use the deterministic equations from (N2) as additional parts of the numerical continuation problem. For example, one may use the numerical methods in (N1) also for the ODE derived in (N2).

There are several advantages of the approach (N1)-(N3). The main advantage is that one works directly with the dynamical objects of interest. For example, if one is interested in the local stochastic fluctuations on a certain time scale depending upon parameters, and one finds an algebraic equation to describe this behaviour in step (N2), then one just focuses the numerical computations upon this task. This generally yields substantial computational- and labour-time savings. Furthermore, the step (N1) can be performed independently, if necessary. This allows for post-processing of models, which have been studied on the deterministic level already via numerical continuation. The disadvantage of (N1)-(N3) is that one actually has to invest a lot of analytical work in (N2) to determine the right equations and in (N3) to determine good discretizations for these equations. However, this effort is mitigated by the fact that this work in (N2)-(N3) has only to be performed once. The situation is similar to the ODE case. One could try to compute bifurcation diagrams for ODEs via MC methods (also called “brute-force” bifurcation diagrams [54]). However, one may also invest quite some effort to find the correct algebraic equations for numerical detection and continuation of bifurcations [53] and then automate the strategy, which leads to very fast and flexible software tools [22, 24].

In this paper, we give a first proof-of-concept that the ideas outlined in (N1)-(N3) have strong potential to work for SPDEs. In particular, we study the Allen-Cahn SPDE (1) with a polynomial nonlinearity


where is the main bifurcation parameter. Next, we briefly summarize the main results and outline the structure of the paper. In Section 2, some background for the deterministic Allen-Cahn PDE is provided. We also motivate the choice (2) of the nonlinearity in more detail. Section 3 very briefly reviews numerical continuation methods for PDEs to fix some notation. In Section 4, we use the toolbox pde2path [74] to compute the main solution branches upon variation of for the Allen-Cahn PDE on an asymmetric rectangular domain. This computation is well-known and yields branch points from a homogeneous solution as well as fold points on the non-trivial branches. In particular, we are interested in those parts of the homogeneous branch and of the first non-trivial branch , which are locally asymptotically stable for the Allen-Cahn PDE.

In Section 5, we discuss the Allen-Cahn SPDE (1). We use a -trace-class Wiener process to define the stochastic forcing term following the standard Hilbert space evolution equation approach to SPDEs. The main goal for the SPDE (1) is to find an efficient numerical method to quantify the stochastic fluctuations depending upon and in phase space regions near deterministic locally asymptotically stable solutions. In this paper, we just focus on the regime of sub-exponential time scales, i.e., before any large deviation effects have occured. To capture the fluctuations, we consider in Section 6 a suitable Lyapunov operator equation, which captures the covariance of the process. Instead of just discretizing the operator equation, we take a more general approach and start from the full Allen-Cahn SPDE in Section 7. The Laplacian is discretized using a finite-difference scheme, the noise term is approximated by using truncation of its spectral representation and the resulting equation is projected onto piecewise constant basis functions. This approach yields a large set of SODEs. We already note that any other spatial discretization technique would also work in the approach (N1)-(N3). In Section 8, we discuss the covariance matrix of a suitable linearized SODE and the associated Lyapunov matrix equation. The covariance matrix for the SODEs is the required approximation for the solution of the full Lyapunov operator equation.

The finite-dimensional Lyapunov equation is solved numerically using iterative solvers in a predictor-corrector framework along parts of the bifurcation branches and in Section 9. The main numerical results are then presented and discussed in Section 9.1 for the homogeneous branch and in Section 9.2 for the non-trivial branch. We summarize the main conclusions in a non-technical way:

  1. The computation of the large covariance matrices can be carried out along a typical bifurcation branch for the Allen-Cahn SPDE within a few minutes on a current standard desktop computer setup. This holds despite the fact that we have not used specialized numerical tools and only used a very simple predictor-corrector scheme. Hence, there is strong potential of the general approach (N1)-(N3) to allow for rapid uncertainty quantification analysis. Furthermore, a comparison between different large-scale linear solvers has been carried out for this problem.

  2. The stochastic fluctuations increase as branch and fold bifurcation points are approached along a stable branch as expected from the theory of critical transitions. The inverse linear scaling for the branch point and the square-root scaling for the fold point were observed in accordance with theoretical predictions.

  3. A numerical comparison between different norms to measure fluctuations, noise truncation levels and various additive as well as multiplicative noise terms has been performed. This shows the flexibility of the approach to characterize the influence of different noise terms on the local dynamics near parameterized steady state branches.

Before proceeding to the main part of the manuscript, we remark that there is a very large number of open problems that could be treated using the philosophy (N1)-(N3) in the context of SPDEs. Here we have just restricted ourselves to a proof-of-concept to present the main ideas clearly. Indeed, there are possible improvements in various different directions, which are currently work in progress. For example, the numerical discretization we have used here is sub-optimal, particularly for pattern-forming problems with spike-layers. One could investigate problems with much more complicated noise and/or aim to answer questions directly arising in larger-scale SPDE models from applications. (N1)-(N3) could also be applied to SPDEs not arising from forcing, but where coefficient functions are random variables. Furthermore, analyzing the stochastic fluctuations numerically very close to bifurcation points requires new ideas to treat the full nonlinear stochastic system in this regime. Overall, the approach to focus directly on dynamical structures for SPDEs in combination with parameter continuation could become another building block in uncertainty quantification.

Acknowledgments: I would like to thank the Austrian Academy of Sciences (ÖAW) for support via an APART fellowship. I also acknowledge support of the European Commission (EC/REA) via a Marie-Curie International Re-integration Grant. Furthermore, I would like to acknowledge the helpful comments of two anonymous referees, which lead to improvements in the presentation of the results.

2 The Cubic-Quintic Allen-Cahn Equation

Let denote a bounded domain. As outlined in the introduction, our basic test problem will be based upon perturbing a version of the Allen-Cahn [2] equation. In general, the parabolic evolution Allen-Cahn PDE for is given by


where is the usual Laplacian, is the potential, is a parameter, indicates the derivative of with respect to the -component and (3) has to be augemented with suitable inital and boundary conditions. We remark that the general Allen-Cahn PDE (3) is a well-studied system for various potentials given by polynomial nonlinearities. There has been a wide range of different analytical and numerical methods proposed to study Allen-Cahn-type PDEs, here we just mention a few examples [1, 4, 65, 76]. For the classical quartic potential leading to a cubic nonlinearity, the Allen-Cahn PDE is also known as the (real) Ginzburg-Landau PDE [36, 5], which is an amplitude equation or normal form for bifurcations from homogeneous states [57, 69]. Furthermore, the Ginzburg-Landau PDE is sometimes referred to as the Nagumo equation [60, 56, 62] developed as a simplified description for waves in nerve axons. Hence, the choice of deterministic Allen-Cahn PDE with a relatively general nonlinearity can definitely be viewed as a standard test problem for new methods. In this paper, we shall restrict to the case


Numerically, the case can also be covered with the techniques developed here, while will probably require more specialized numerical ideas. Basically we focus on here to have a sufficiently challenging numerical setup. The stationary problem associated to (3), for , with Dirichlet boundary conditions, is given by


As a potential we are going to fix the following function


Since the nonlinearities in are of cubic and quintic types, one also refers to (3) or (5) with the choice (6) as the cubic-quintic Allen-Cahn (cqAC) equation. In fact, it is also very common to add a quintic term to the Ginzburg-Landau case, see for example [21, 42, 75], which again justifies our choice of the qcAC equation as a basic test example.

3 Numerical Continuation Basics

It is natural to ask, which solutions the stationary cubic-quintic Allen-Cahn equation (5) has under the variation of . Numerical continuation methods are a natural tool to address this problem. First, one uses a spatial discretization, e.g., finite elements or finite differences to convert (5) into a nonlinear algebraic system of equations


for some , where is usually very large for PDE problems. The idea of numerical continuation is to trace out a curve of solutions where parametrizes the the curve, for example by arclength. There are theoretically well-established predictor-corrector schemes available [34, 70] to compute the curve . The main idea is that knowing a solution point then one can compute with , for sufficiently small, very efficiently. Indeed, first one uses a prediction step, for example by tangent prediction


where is the tangent vector to the curve at . Then ones uses a correction step, i.e., one has to solve


iteratively with starting initial guess . If the initial guess is sufficiently close to a true solution, which can be guaranteed under suitable non-degeneracy hypotheses, then an iterative solver converges very quickly to . At degenerate bifurcation points, one may use branch switching techniques to track the bifurcating curves beyond this point [53].

4 Deterministic Continuation for the cqAC Equation

Figure 1 shows the results of numerical continuation for the cqAC equation (5). The results were obtained from a standard example in the continuation software pde2path [74] on a rectangle


The continuation was started with the homogeneous solution , which exists for all . Increasing from the starting value , we focused on the first three bifurcation points. Those points are approximately obtained at


We remark that it is known for the square domain that the bifurcation points and coincide, which is one reason to consider the rectangle to split up this numerical difficulty.

Figure 1: (a) Bifurcation diagram of the cqAC equation (5)-(6) obtained from four continuation runs in pde2path; the vertical axis shows the -norm of the solution. The thick curves (black and blue) are locally stable, i.e., for the linearized problem the spectrum is contained in the left-half of the complex plane. The three simple branch/bifurcation points (11) are marked with circles. There are also three fold points (12) on the non-trivial bifurcation branches (blue/red/magenta). The five formally defined parameter regimes (13) are delineated by vertical dashed green lines. The three black dots indicate to which points the solutions in diagrams (b)-(d) are associated. (b)-(d) Solution plots for parameter values on the upper parts of the non-trivial branches.

Up to , the homogeneous branch is locally asymptotically stable for the time-dependent cqAC equation (3). The linearized system along the homogeneous branch has at least one eigenvalue with positive real part for and only eigenvalues with negative real part for . Branch switching at , , in pde2path leads to three subcritical bifurcation branches , which each re-acquire one negative eigenvalue at further fold points located at


In particular, the branch restabilizes and the system is bistable for . We note the numerical results clearly indicate that is globally stable for , while the top part of the branch is the only attractor for .

Remark: At the bifurcation points, symmetries of the cqAc equation considered here ( or or ) imply the existence of additional bifurcation branches with symmetric solutions not shown in Figure 1. In particular, throughout this manuscript we only select and analyze numerically one of the symmetric branches. For more details on symmetry-breaking and using numerical continuation to find starting solutions without homotopy continuation see, e.g., [51].

We partition the parameter space into five disjoint intervals, which we are going to focus on throughout this manuscript


Basically, the regime is globally monostable, is bistable, has one stable and one simple homogeneous unstable branch, while and consider the interaction between one and two non-trivial unstable branches and , respectively, with the stable branch and the homogeneous branch .

5 The Stochastic Allen-Cahn Equation

We return to the general Allen-Cahn equation (3). One natural approach to include noise is to replace the PDE, on a very formal level, by a stochastic partial differential equation (SPDE) given by


where is a space-time-dependent stochastic process and is a suitable mapping. The Allen-Cahn SPDE has gained - unsurpringly considering the discussion in Section 2 - considerable attention in various analytical studies, see e.g. [26, 12, 16, 38, 45, 77]. Furthermore, numerical methods for SPDEs have frequently been tested on 2D Allen-Cahn-type SPDEs and related equations [19, 46, 55, 66, 71]. Hence, (14) can also be viewed as a standard test problem for the SPDE case.

To make the formal SPDE (14) precise, we shall adopt the framework of mild solutions considerd in [64]. Consider a probability space and a Hilbert space ; for practical purposes we may always think of or here. Denote by the inner product in and by the associated norm. Let be a symmetric non-negative linear operator and suppose is of trace-class . Then there exists a complete orthonormal system for the Hilbert space such that

where is a non-negative bounded summable of eigenvalues. Then one may define a -Wiener process by


where are independent one-dimensional Brownian motions. The series (15) converges in ; see also [64, p.86-89]. One finds that and for all we have so that can be viewed as the covariance operator. Furthermore, define and let , where denotes the space of Hilbert-Schmidt operators. In this case one may consider the stochastic Allen-Cahn equation for as an -valued stochastic evolution equation


which can interpreted in the usual precise integral form [64] by stochastic integrals; note that in (16), we understand in the Nemytskii operator sense as discussed in [64, Sec.7.2] and mild solutions exist under suitable local Lipschitz conditions on . Furthermore, we shall always assume that the initial condition is deterministic and lies in a suitable subspace of , e.g.  as discussed in [64].

The solution theory for general SPDEs is quite involved and is still being developed for certain cases. However, we are going to work with a linearized version of (16) for which a suitable theory of mild solutions exists; see Section 6.

6 The Abstract Lyapunov Equation

We have seen in Section 4 that the deterministic cqAC PDE has several locally asymptotically stable stationary solutions. It is then a natural question to ask, what influence different forms of the noise term have on the fluctuations around these stationary solutions. Here we are interested in fluctuations on a subexponential time-scale, i.e., before the regime of large deviation theory [27] is reached. This implies that we restrict to either

  • sufficiently small norm of , or

  • sufficiently small eigenvalues , or

  • sufficiently small time final such that no large deviations have occurred.

In this regard, it is important to observe that (R1)-(R3) are related. For example, may absorb a constant positive multiplicative prefactor in either into by defnining , or into the eigenvalues associated to by defining a new operator . Furthermore, one may apply a change of time scale in (16) via for a new time and a positive scale factor to relate the condition (R3) to (R1)-(R2).

In one (or more) of the cases (R1)-(R3), it natural to consider a linearized form of (16) around a deterministic stationary solution for a parameter value . The linear SPDE is given by


where , the maps are linear operators given by


and denotes the second derivative of with respect to the -component, i.e., we have


where just denotes the identity operator. Note carefully that the spectrum of , denoted here by , indicates whether is a locally stable solution of the deterministic cqAC PDE (3). If is properly contained in the left-half complex plane, then we have local asymptotic stability. This is the main case we are going to focus on here; recall that the locally stable stationary solutions of the cqAC PDE have been computed solving (5) by numerical continuation in Section 4 and are shown as the thick curves in Figure 1.

The linear SPDE (17) is the direct analog of the stochastic ordinary differential equation for an Ornstein-Uhlenbeck (OU) process. As a solution concept, we consider mild solutions [64, Ch.7] given by


where is the semigroup generated by . It is helpful to denote the stochastic integral in (20) as

is also called the stochastic convolution, which has the associated covariance operator [64, Thm. 5.2]


where denotes the adjoint operator of and denotes the adjoint semigroup of . Computing provides the main relevant information on the local fluctuations of the OU-process. Since the OU process differs only in the quadratic and higher-order terms from the local dynamics around the stationary solution of the cqAC SPDE, one may show that on subexponential time scales, the operator provides the leading-order fluctuations around ; see also Appendix A.2. Using (15), one can write the OU process solution also as


where the detailed definition of the integral is given in [63, Sec. 2.2]. Using a standard result [63, Prop. 2.2], it is not difficult to see that the series in (22) is convergent in . Furthermore, is a Gaussian random variable with mean and covariance operator . Indeed, if we are in the case where is locally stable for the deterministic cqAC PDE then we have

for some constant and , where . Therefore, as exponentially fast, so we may formally neglect the first term . In this case, we simply have and there exists a unique invariant Gaussian measure with mean zero and covariance operator [63, Thm. 2.34]; again the convergence to is expected to be exponentially fast so we may negelect the initial fast transients and just focus on . is a linear continuous operator that satisfies for all and . Furthermore, it is crucial to note that satisfies a Lyapunov equation [63, Lem. 2.45] given by


for any . Observe that (23) has the same algebraic structure of the classical Lyapunov equation associated with linear stochastic ordinary differential equations (SODEs) [11, Ch.5]. The result (23) can also be extended to apply in Banach spaces [32, Sec.4].

7 A Finite-Difference Galerkin Scheme

One possibility to compute the operator is to start directly from the abstract Lyapunov equation (23) and derive suitable numerical methods for it. However, this studies (23) in isolation from the other numerical algorithms one may develop for the deterministic and stochastic cqAC equation. Here we aim to determine a finite-dimensional approximation to by starting from the SPDE (16).

As a first step, we use a spatial finite-difference discretization for the SPDE (16) based upon the one-dimensional scheme presented in [68], in combination with a spectral approximation of the noise [71]. The reason for this combination will be explained below, after we have introduced the scheme. Let for and denote the points in a regular spatial mesh of edge lengths for the rectangle , i.e., we take


which yields points and interior points. Denote the approximation of on the mesh by . We approximate the Laplacian on the interior points, assuming Dirichlet boundary conditions, by the usual five-point stencil given by

The nonlinear terms are approximated by evaluation at the mesh points


where the -dependence of the deterministic nonlinear part will be omitted in the notation for simplicity. Next, denote by


the -normalized indicator function of the rectangle . To approximate the covariance operator consider, the expansion


where we may select to control the approximation accuracy of the noise term. Then we consider the inner product of the cqAC SPDE with the indicator function (26) to obtain


Now one may consider each of the different terms for and note that


Similarly, we have


For the noise term, we use a truncation level as in (27), and obtain


where are scalars which can be pre-computed in an offline step. We use the notation , where is an index labelling the interior vertices of the mesh and which ranges between and , consider the vector and define


so that is just a vector-valued map and is a matrix-valued map. Hence, we may approximate the SPDE with a -dimensional system of SODEs given by


where is a vector of independent identically distributed (iid) Brownian motions. Note that we have now two different ways to control the approximation level. The spatial mesh can be refined by increasing and/or while the spectral approximation of the noise can be improved by increasing . Of course, many other spatial discretization schemes would lead to a system of SODEs. The next step starts from the level of the SODEs (33) to approximate the covariance operator of the problem.

8 Lyapunov Equations for SODEs

We start from the SODE (33), where we are interested primarily in the case where the noise level is suitably bounded, i.e., for all in some compact set of phase space and is a fixed constant. This is quite a natural viewpoint as very large noise terms may indicate that the model might have not been derived to suitable physical modelling accuracy in the first place. Here we briefly recall the idea to approximate certain covariance operators associated to (33) along bifurcation curves developed in [48].

Suppose that the associated ODE for , given by


has a hyperbolic stable steady state for a given range of parameter values. For the cqAC SPDE case we consider here, we may think of the discretized versions of the stable parts of the branches and computed in Section 4. For example, for the homogeneous branch in the parameter range we have the stable solution for all . If we suppose that and Taylor expand (33) to lowest order at we obtain


where is the usual Jacobian matrix and is the lowest order term of the noise. Equation (35) is a -dimensional Ornstein-Uhlenbeck (OU) process [30] with covariance matrix

where is the fundamental solution of . Differentiation shows that satisfies the ODE


Since is a hyperbolic stable equilibrium point, it follows [8] that the eigenvalues of the linear operator are , where are the eigenvalues of . Therefore, the covariance ODE (36) has a globally stable equilibrium solution defined by


which is obtained by solving


Note carefully that (38) is a Lyapunov equation and is a finite-dimensional approximation to the covariance operator discussed in Section 6. We denote the solution of (38) sometimes also as to stress the dependence upon the parameter and the associated state .

It is standard [11] to define a covariance neighbourhood


where is a parameter that can be interpreted as a probabilistic confidence level. Essentially, yields the information about local fluctuations and important noise directions near . Hence, we would like to solve for along the entire branch of steady states obtained via continuation . The algebraic equation (38) is a uniquely solvable Lyapunov matrix equation [28]. If the dimension of the problem is moderate, direct numerical methods for (38) work quite well. A well-known direct solution algorithm is the Bartels-Stewart algorithm [7, 33].

However, we remark that for our case new aspects arise since we want to solve (38) along an entire steady state branch . Indeed, most numerical continuation algorithms require an approximation of the Jacobian matrix to compute a point starting from . Therefore, the matrix is available at each continuation step. Furthermore, assembling the matrix at a given point is relatively cheap numerically if is of moderate size. Solving (38) at gives a matrix . If is small then , or a tangent prediction starting from , is already an excellent initial guess to find . Hence, except for the first point on the equilibrium curve, we always have an initial guess available for an iterative method of solving the Lyapunov equation. This was exploited in [48] via a Gauss-Seidel scheme by re-writing the Lyapunov equation as a linear system.

9 Continuation of the Lyapunov Equation for cqAC

In this paper, we also re-write the Lyapunov equation (38) as a large (sparse) linear system using the Kronecker product in the usual way


where Id is the identity matrix of size and converts a matrix to a vector by stacking it column-by-column, e.g., . Here we are going to use standard iterative solvers for sparse linear systems to continue the Lyapunov equation as implemented in MatLab. In particular, we use bicgstab (a stabilized biconjugate gradient method, our default choice here), gmres (a generalized minimum residual method) and qmr (a quasi-minimal residual method). We use as an initial guess the previous solution of the Lyapunov equation on the continuation branch; if we have no solution available, the zero vector is used as an initial guess.

To obtain (40), we also have to make a choice of the mesh parameters described in Section 7. Note that we work on a regular grid on the rectanglular domain


if and only if , which implies upon fixing that . Hence, to obtain the case , we must require that is a natural number. One good choice is and so , which will be our standard choice here. Recall again that here we do not aim for optimal numerical performance of the algorithm. Indeed, it is important to note that the method we use here to continue the Lyapunov equation is numerically sub-optimal in several different ways:

  1. There are more efficient methods for iterative solution of matrix Lyapunov equations available that take direct advantage of the special matrix-equation structure, e.g. ADI-type methods would be an option [9, 28].

  2. Instead of using the previous solution as an initial guess, a full numerical continuation scheme may also employ other predictors, e.g. a tangent predictor [34].

  3. Using a regular mesh is likely to be sub-optimal in many cases. One could use a suitable finite-element method spatial discretization together with adaptive meshing for the spatial discretization, which is going to improve the performance in many situations.

If, despite using a sub-optimal numerical algorithm as explained in (S1)-(S3), we can still compute the Lyapunov equation efficiently on a standard desktop hardware platform111Basic details of the desktop computer setup used: Intel Core i5-4430 CPU @ 3.00GHz processor (quadcore), Kingston 16GB system memory, WDC WD10EZRX-00L 1TB harddrive, Ubuntu 12.04.4 LTS operating system, MatLab R2013a computational platform. with our basic continuation setup, then we have achieved a valid proof-of-concept. In fact, this leads to the conjecture that the computation time can be reduced potentially by one or more orders of magnitude if all available ideas in (S1)-(S3) are explored.

We use the domain and set as our Hilbert space. Furthermore, as a basic test noise we consider a decaying noise term with eigenvalues for given by


where is a parameter, which can be used to tune the noise level and is a vector with increasing non-negative entries to be specified. As a basis , we fix


As default parameters for the continuation of the Lyapunov equations we set and consider an additive noise as a first test case so that . Regarding the linear system solvers, we fix the maximum iteration number for the iterative linear system solvers to and the tolerance to .

9.1 The Lower Branch

For the first numerical test case, we focus on the homogeneous stationary solution branch , which is locally asymptotically stable for with for the deterministic system. Three different noise terms were considered with eigenvalues given by (42) associated to three different truncation levels , more precisely, we have taken for .

Figure 2: The computation time and number of iterations for the Lyapunov operator computation for the lower homogeneous branch with are shown. Each circle represents a point on the continuation curve, where (38) has been solved using the bicgstab-method. The linear interpolation between the points is for visualization purposes. (a) Parameter plotted against computation time (in seconds). The cyan (), green () and black () are different truncation levels for the noise. (b) Number of iterations required for bicgstab to reach a tolerance of less than . The color coding is as for (a).
Figure 3: (a) All norms in this figure are maximum norms . The deterministic homogeneous steady state solution (thick black curve) is shown, which also corresponds to the zero steady state of the discretized system. Furthermore, the norm of the covarince matrix is plotted to indicate the size of the fluctuations as measured by the linearized system; there are actually three different neighbourhoods shown in cyan (), green () and black () corresponding to different truncation levels for the noise. The inset in (a) shows a semi-logarithmic plot of the upper curves (thin curves) in comparison to the scaling law (dashed black thick curve). (b) Zoom of part (a) near the upper part of the covariance neighbourhood to show that the curves given by for are very close. (c) Eigenvalues used for the noise, i.e., for there are two eigenvalues (in cyan), for there are fo and so ur eigenvalues (the cyan ones and the green ones), while for also the black ones are taken into account.

The results of the computations are shown in Figures 2-3. Figure 2 considers the basic computational performance for the large sparse linear systems solution to find . The computational time and the number of iterations are very large in the first step, as we do not have a good initial guess available and just start with the zero vector. Afterwards, each continuation point takes approximately between 20 and 30 seconds on a standard desktop computer so the entire branch of 17 points is processed within a few minutes. We conjecture that taking (S1)-(S3) into account, the computation time can be reduced to a few seconds for the entire branch on the same hardware. In any case, even for the sub-optimal computational framework, we can quickly analyze the local influence of noise and quantify the uncertainty of the local system dynamics for different parameter values.

Figure 4: Comparison of different norms for the vector of variances , which is obtained as the diagonal of . As usual, the horizontal axis in (a)-(c) shows the parameter values for , where the homogeneous zero branch (thick black) is locally asymptotically stable. The figures (a)-(c), respectively crosses, dots and circles, correspond to , and .

and so Figure 3(a) shows the influence of noise in the maximum norm. Note that the solution corresponds to the finite-dimensional approximation , which is a steady state of the ODE (34). The results show that the stochastic neighbourhood defined by the maximum norm of the covariance matrix widens as the bifurcation point is approached. This is due to critical slowing down, which is an effect well-studied in certain regimes near bifurcation points for SODEs [49]. In particular, the increase in the (co-)variance can be used as an early-warning sign for the upcoming bifurcation point. For SPDEs, the scaling of the covariance operator has been investigated recently using analytical techniques in one-dimensional normal form-type SPDEs for pattern formation near branch points [35], where it is shown that we formally expect a scaling


in certain modes and , where is a constant. Since we are taking the maximum norm here, the result (44) does not apply directly but we still see an increase in that shows a qualitatively similar in the inset in Figure 3(a), where we plot (44) with . However, there are differences near the bifurcation point, which still have to be explored in future work; see also [35] where this issue has been partially discussed. Of course, very close to a bifurcation point, one has to take into account the full nonlinear dynamics of the problem and the linearization approach via Ornstein-Uhlenbeck processes is no longer valid; we also leave this as an extension for future work.

Figure 5: Visualization of the stochastic neighbourhood calculated from the covariance matrix for with for . Each point in the grid corresponds to a point and on the vertical axis labelled , we show the stationary solution (a green plane) as well as the associated variances for each point (the surface with coloring green/yellow/red) and (the surface with coloring green/cyan/blue). The mesh is overlayed on all three objects. Note that it is difficult to visualize the full evolution of the associated variances when varying the parameter in a three-dimensional plot so this paper contains as a supplementary material a movie, which displays this figure for .

Figure 3(b) shows a zoom near the covariance neighbourhood boundary curves defined by to illustrate that there is only a very small difference between the three noise truncation levels . In fact, the differences are very close to the numerical precision used for the linear solver . This result is expected since the largest stochastic fluctuation, when measured in the , is dominated by the mode with the largest eigenvalue. The largest eigenvalue is included in all three cases as shown in Figure 3(c).

Figure 4 compares the growth of the stochastic fluctuations in different norms. In this case, we focused only on the vector of variances


Note carefully that the scale difference in Figure 4 is over three orders of magnitude, so one has to be very careful to define a threshold for covariance neighbourhoods only when it is clear, which norm is used to evaluate this threshold. The most important aspect of Figure 4 is that in all the different norms, an early-warning sign, based upon the scaling law of the covariance neighbourhoods before the bifurcation point, is clearly visible.

Figure 6: Computation time for Lyapunov operator computation for the lower homogeneous branch with for three different iterative linear solvers. Each circle represents a point on the continuation curve, where (38) has been solved using bicgstab (cyan), gmres (green/magenta) and qmr (black). For gmres, there is also the option to restart after a fixed number of inner iterations, which has been set to 10 (magenta) and 0 (green) respectively. The linear interpolation between the points is for visualization purposes. The parameter is plotted against computation time (in seconds) in (a) and (b). Note that the gmres method in (b) has a lot longer computation time.

Of course, one may also consider a visualization of the results in three dimensions. In Figure 5 the stationary solution is shown together with a stochastic neighbourhood defined via the variances . More precisely, each grid point has an associated solution value and an associated variance ; of course, we plot also . One clearly sees that the fluctuations have been computed for the Dirichlet Laplacian with zero boundary conditions so that the fluctuations are largest at the center of the domain. Figure 5 shows the case when ; in the supplementary material a movie is provided showing the growth of the stochastic neighbourhood in three dimensions upon variation of ; for a brief comparison to Monte-Carlo simulations of the full nonlinear system see Appendix A.1.

Figure 7: Computation of the upper locally asymptotically stable branch for . (a) All norms in this figure are maximum norms . We do not show the initial part in region near the fold point in this part of the figure; see Figure 8. The deterministic steady state solution on the branch (thick blue curve) from the deterministic continuation run in Figure 1 is shown together with the unstable branch (red). The branches and correspond to steady states of the discretized system (34). Furthermore, the norm of the covarince matrix is plotted to indicate the size of the fluctuations as measured by the linearized system (thin blue curves). There are actually four different neighbourhoods shown for additive noise, eigenvalues (42) for . (b) Computation time along the entire upper branch. (c) Number of iterations . In (b)-(c) each dot represents a point on the continuation curve, where (38) has been solved using the bicgstab-method. Only the computation for is shown, the other three cases yield very similar results.

The performance of different standard iterative linear solvers for the large sparse linear system (40) has also been investigated. The lower homogeneous branch has again been used as a test case. We used bicgstab, gmres and qmr in their MatLab implementation222Version: MatLab R2013a with tolerance . As before, we use the eigenvalues (42) now with for . Figure 6 shows a comparison in the computation time for each point on the bifurcation curve. The results show that there are certainly differences in the computation time for the different methods and that the bicgstab method is the fastest in our context while gmres is the slowest for this problem. Given the different performances observed here for various iterative linear system solvers for the Lyapunov equations we have to solve, we can conjecture that using a specialized algorithm for Lyapunov equations may generate substantial additional speed-up as indicated in (S1) above.

9.2 The Upper Branch

In the next numerical continuation run, we focus on the upper locally asymptotically stable steady state branch computed in Section 4; see Figure 1. The algorithmic parameters for the linear system solvers are kept as for the homogeneous branch computations from Figures 2-5 now with for and eigenvalues given by (42) for additive noise. Furthermore, we recall that there is a locally asymptotically stable non-trivial solution branch in the parameter regime as shown in Figure 1. The branch changes stability at the fold point .

As before, we use the eigenvalues (42) now with for .

Figure 8: Computation of the upper locally asymptotically stable branch for , now with a focus on the region near the fold point at . (a) All norms in this figure are maximum norms . The color coding of (a) is the same as in Figure 7(a). There are four different neighbourhoods shown for additive noise, eigenvalues (42) for . (b) Computation of the scaling laws for the growth of the norm as the bifurcation point is approached. The four thin blue curves are just a double-logarithmic plot of the curves from (a) with the bifurcation point shifted to zero. The thick black dashed line has slope , while the thick black solid line has slope ; the latter scaling fits the observed scaling law very well.

Figure 7(a) shows the continuation run with the branch and a focus on , i.e., bounded away from the fold point. Four stochastic neighbourhoods computed for different noise strengths are shown. If we increase , then the deterministic stability increases. Since the additive noise is kept constant, this explains the shrinking of the four neighbourhoods. Note carefully, that the neighbourhoods will interact with the unstable branch . The precise interaction threshold depends upon the nonlinearity but we can observe that there is a competition effect: the norm of the solution on the unstable branch increases while the size of the stochastic neighbourhoods decreases. In the current setup, we can observe that the noise-induced effect is going to dominate eventually for sufficiently large and we expect frequent noise-induced excursions on a sub-exponential time-scales across the unstable branch for individual sample paths.

Figures 7(b)-(c) show the computational performance along the branch . As we start with the zero vector as an initial guess for the first value smallest value of near , the iterative method is initially slow but again increases performance after a few points on the continuation branch as in Figure 2. We also observe from the results that it seems to be computationally faster to compute the Lyapunov equation in a regime, where the leading eigenvalue of the linearized deterministic system is smaller, i.e., in a more stable regime.

Remark: In Figures 7-8 we have increased the noise level by several orders of magnitude to visualize the neighbourhoods better in comparison to the order one scale of the deterministic bifurcation branch. Note that the Lyapunov equation is invariant under a change of scale, i.e., a multiplicative factor for in front of the noise term defined by enters as a scale factor of for the covariance in; see the Lyapunov equation (40). This avoids having to add several zoom-ins near the bifurcation branch in Figures 7-8.

Figure 8 focuses on the fold point. In Figure 8(a), we consider and illustrate the growth of the stochastic fluctuations near the bifurcation point. Again, we are interested in scaling laws. For example, the goal is to determine the exponent in


Figure 8(b) shows that does not yield a good fit near the fold in contrast to the case near the bifurcation from the homogeneous branch at . The fit for is a lot better and this corresponds precisely to the theory expected from SODEs. Indeed, near transcritical and pitchfork bifurcations one finds for SODEs , while for the fold one obtains [49]. Hence, the result in Figure 8(b) is a natural generalization from the SODEs results as expected from the recent SPDEs results on early-warning signs in [35].

Figure 9: Results for different types of scalar-valued noise terms for fluctuations around the upper locally asymptotically stable deterministic bifurcation branch with . All norms in this figure are maximum norms (a) The deterministic homogeneous steady state solution on the branch (thick blue curve) from the deterministic continuation run in Figure 1 is shown together with the unstable branch (red). Three different noise terms are considered and the associated neighbourhoods are plotted: additive noise (thin blue curves), norm-dependent noise scaling (47) (thin cyan curves) and norm-dependent noise amplification/reduction (48) (thin magenta curves). The value of has been fixed and the eigenvalues (42) with truncation have been used. (b)-(c) Zooms of different parts of the curves in (a).

In the next numerical continuation run we investigate different noises. The most natural first generalization from purely additive noise is to consider scalar-valued noise terms , where does depend upon , i.e., relatively simple multiplicative noise terms. The choices we make here are purely for test purposes and many other situations may arise in practical applications. We consider additive noise and the following two cases


The term (47) models a scaling effect of the noise depeding upon the current maximal value of the solution. The term (48) is a shifted multiplicative noise term. The results in Figure 9 show that multiplicative noise can indeed crucially change the relative behaviour in different regimes. For example, if so that we are near the fold point, then the stochastic neighbourhood generated by the noise-amplification term (47) is smaller than the additive noise neighbourhood, while the situation is interchanged for , i.e., there is a change in the intermediate regimes when passes through . Furthermore, the results show that it is a-priori not easy to see from the structure of the noise term, which stochastic neighbourhood is larger when two different types of multiplicative noise are compared. Hence, the method proposed here can automate this process, even as a fast post-processing step for deterministic bif