Numerical Continuation and SPDE Stability
for the 2D CubicQuintic AllenCahn Equation
Abstract
We study the AllenCahn equation with a cubicquintic nonlinear term and a traceclass 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 finitedifferences and spectral noise approximation obtaining a finitedimensional 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 suboptimal computational setup, we can quantify the subexponentialtimescale 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 spatiotemporal stochastic systems.
Keywords: AllenCahn, SPDE, cubicquintic GinzburgLandau, Lyapunov equation, numerical
continuation, spectral noise approximation, uncertainty quantification, finite differences,
iterative linear solvers, critical transitions, bifurcation diagram.
1 Introduction
Spatiotemporal 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 spatiotemporal 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 AllenCahn SPDE, which can formally be written as
(1) 
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 AllenCahntype SPDE are discussed below.
Currently, a standard approach to study the dynamics of (1) numerically is centered around MonteCarlo (MC) algorithms and their variations. The strategy is to use forward integration of the SPDE. This requires three different components:

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

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

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 quasiMonteCarlo methods [18, 61] based on discrepancy theory or multilevelMonteCarlo [39, 31] based upon a coarse/finesampling hierarchy. However, there are some drawbacks when following an MCbased 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 laborintensive. 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.12].
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 predictorcorrector 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:

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.

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 prefactor and exponent in Kramers’ law governing metastable gradient systems.

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 labourtime savings. Furthermore, the step (N1) can be performed independently, if necessary. This allows for postprocessing 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 “bruteforce” 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 proofofconcept that the ideas outlined in (N1)(N3) have strong potential to work for SPDEs. In particular, we study the AllenCahn SPDE (1) with a polynomial nonlinearity
(2) 
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 AllenCahn 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 AllenCahn PDE on an asymmetric rectangular domain. This computation is wellknown and yields branch points from a homogeneous solution as well as fold points on the nontrivial branches. In particular, we are interested in those parts of the homogeneous branch and of the first nontrivial branch , which are locally asymptotically stable for the AllenCahn PDE.
In Section 5, we discuss the AllenCahn SPDE (1). We use a traceclass 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 subexponential 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 AllenCahn SPDE in Section 7. The Laplacian is discretized using a finitedifference 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 finitedimensional Lyapunov equation is solved numerically using iterative solvers in a predictorcorrector 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 nontrivial branch. We summarize the main conclusions in a nontechnical way:

The computation of the large covariance matrices can be carried out along a typical bifurcation branch for the AllenCahn 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 predictorcorrector scheme. Hence, there is strong potential of the general approach (N1)(N3) to allow for rapid uncertainty quantification analysis. Furthermore, a comparison between different largescale linear solvers has been carried out for this problem.

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 squareroot scaling for the fold point were observed in accordance with theoretical predictions.

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 proofofconcept 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 suboptimal, particularly for patternforming problems with spikelayers. One could investigate problems with much more complicated noise and/or aim to answer questions directly arising in largerscale 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 MarieCurie International Reintegration 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 CubicQuintic AllenCahn Equation
Let denote a bounded domain. As outlined in the introduction, our basic test problem will be based upon perturbing a version of the AllenCahn [2] equation. In general, the parabolic evolution AllenCahn PDE for is given by
(3) 
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 AllenCahn PDE (3) is a wellstudied system for various potentials given by polynomial nonlinearities. There has been a wide range of different analytical and numerical methods proposed to study AllenCahntype PDEs, here we just mention a few examples [1, 4, 65, 76]. For the classical quartic potential leading to a cubic nonlinearity, the AllenCahn PDE is also known as the (real) GinzburgLandau PDE [36, 5], which is an amplitude equation or normal form for bifurcations from homogeneous states [57, 69]. Furthermore, the GinzburgLandau 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 AllenCahn 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
(4) 
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
(5) 
As a potential we are going to fix the following function
(6) 
Since the nonlinearities in are of cubic and quintic types, one also refers to (3) or (5) with the choice (6) as the cubicquintic AllenCahn (cqAC) equation. In fact, it is also very common to add a quintic term to the GinzburgLandau 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 cubicquintic AllenCahn 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
(7) 
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 wellestablished predictorcorrector 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
(8) 
where is the tangent vector to the curve at . Then ones uses a correction step, i.e., one has to solve
(9) 
iteratively with starting initial guess . If the initial guess is sufficiently close to a true solution, which can be guaranteed under suitable nondegeneracy 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
(10) 
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
(11) 
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.
Up to , the homogeneous branch is locally asymptotically stable for the timedependent 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 reacquire one negative eigenvalue at further fold points located at
(12) 
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 symmetrybreaking 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
(13) 
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 nontrivial unstable branches and , respectively, with the stable branch and the homogeneous branch .
5 The Stochastic AllenCahn Equation
We return to the general AllenCahn 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
(14) 
where is a spacetimedependent stochastic process and is a suitable mapping. The AllenCahn 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 AllenCahntype 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 nonnegative linear operator and suppose is of traceclass . Then there exists a complete orthonormal system for the Hilbert space such that
where is a nonnegative bounded summable of eigenvalues. Then one may define a Wiener process by
(15) 
where are independent onedimensional Brownian motions. The series (15) converges in ; see also [64, p.8689]. 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 HilbertSchmidt operators. In this case one may consider the stochastic AllenCahn equation for as an valued stochastic evolution equation
(16) 
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].
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 timescale, 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
(17) 
where , the maps are linear operators given by
(18) 
and denotes the second derivative of with respect to the component, i.e., we have
(19) 
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 lefthalf 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 OrnsteinUhlenbeck (OU) process. As a solution concept, we consider mild solutions [64, Ch.7] given by
(20) 
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]
(21) 
where denotes the adjoint operator of and denotes the adjoint semigroup of . Computing provides the main relevant information on the local fluctuations of the OUprocess. Since the OU process differs only in the quadratic and higherorder 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 leadingorder fluctuations around ; see also Appendix A.2. Using (15), one can write the OU process solution also as
(22) 
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
(23) 
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 FiniteDifference 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 finitedimensional approximation to by starting from the SPDE (16).
As a first step, we use a spatial finitedifference discretization for the SPDE (16) based upon the onedimensional 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
(24) 
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 fivepoint stencil given by
The nonlinear terms are approximated by evaluation at the mesh points
(25) 
where the dependence of the deterministic nonlinear part will be omitted in the notation for simplicity. Next, denote by
(26) 
the normalized indicator function of the rectangle . To approximate the covariance operator consider, the expansion
(27) 
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
(28) 
Now one may consider each of the different terms for and note that
(29) 
Similarly, we have
(30) 
For the noise term, we use a truncation level as in (27), and obtain
(31) 
where are scalars which can be precomputed 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
(32) 
so that is just a vectorvalued map and is a matrixvalued map. Hence, we may approximate the SPDE with a dimensional system of SODEs given by
(33) 
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
(34) 
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
(35) 
where is the usual Jacobian matrix and is the lowest order term of the noise. Equation (35) is a dimensional OrnsteinUhlenbeck (OU) process [30] with covariance matrix
where is the fundamental solution of . Differentiation shows that satisfies the ODE
(36) 
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
(37) 
which is obtained by solving
(38) 
Note carefully that (38) is a Lyapunov equation and is a finitedimensional 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
(39) 
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 wellknown direct solution algorithm is the BartelsStewart 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 GaussSeidel scheme by rewriting the Lyapunov equation as a linear system.
9 Continuation of the Lyapunov Equation for cqAC
In this paper, we also rewrite the Lyapunov equation (38) as a large (sparse) linear system using the Kronecker product in the usual way
(40) 
where Id is the identity matrix of size and converts a matrix to a vector by stacking it columnbycolumn, 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 quasiminimal 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
(41) 
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 suboptimal in several different ways:

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].

Using a regular mesh is likely to be suboptimal in many cases. One could use a suitable finiteelement 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 suboptimal numerical algorithm as explained in (S1)(S3), we can still compute the Lyapunov equation efficiently on a standard desktop hardware platform^{1}^{1}1Basic details of the desktop computer setup used: Intel Core i54430 CPU @ 3.00GHz processor (quadcore), Kingston 16GB system memory, WDC WD10EZRX00L 1TB harddrive, Ubuntu 12.04.4 LTS operating system, MatLab R2013a computational platform. with our basic continuation setup, then we have achieved a valid proofofconcept. 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
(42) 
where is a parameter, which can be used to tune the noise level and is a vector with increasing nonnegative entries to be specified. As a basis , we fix
(43) 
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 .
The results of the computations are shown in Figures 23. 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 suboptimal computational framework, we can quickly analyze the local influence of noise and quantify the uncertainty of the local system dynamics for different parameter values.
and so Figure 3(a) shows the influence of noise in the maximum norm. Note that the solution corresponds to the finitedimensional 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 wellstudied in certain regimes near bifurcation points for SODEs [49]. In particular, the increase in the (co)variance can be used as an earlywarning sign for the upcoming bifurcation point. For SPDEs, the scaling of the covariance operator has been investigated recently using analytical techniques in onedimensional normal formtype SPDEs for pattern formation near branch points [35], where it is shown that we formally expect a scaling
(44) 
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 OrnsteinUhlenbeck processes is no longer valid; we also leave this as an extension for future work.
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
(45) 
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 earlywarning sign, based upon the scaling law of the covariance neighbourhoods before the bifurcation point, is clearly visible.
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 MonteCarlo simulations of the full nonlinear system see Appendix A.1.
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 implementation^{2}^{2}2Version: 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 speedup 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 25 now with for and eigenvalues given by (42) for additive noise. Furthermore, we recall that there is a locally asymptotically stable nontrivial solution branch in the parameter regime as shown in Figure 1. The branch changes stability at the fold point .
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 noiseinduced effect is going to dominate eventually for sufficiently large and we expect frequent noiseinduced excursions on a subexponential timescales 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 78 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 zoomins near the bifurcation branch in Figures 78.
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
(46) 
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 earlywarning signs in [35].
In the next numerical continuation run we investigate different noises. The most natural first generalization from purely additive noise is to consider scalarvalued 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
(47)  
(48) 
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 noiseamplification 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 apriori 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 postprocessing step for deterministic bifurcation d