A multilevel Monte Carlo finite element method for the stochastic Cahn–Hilliard–Cook equation

In this paper, we employ the multilevel Monte Carlo finite element method to solve the stochastic Cahn–Hilliard–Cook equation. The Ciarlet–Raviart mixed finite element method is applied to solve the fourth-order equation. In order to estimate the mild solution, we use finite elements for space discretization and the semi-implicit Euler–Maruyama method in time. For the stochastic scheme, we use the multilevel method to decrease the computational cost (compared to the Monte Carlo method). We implement the method to solve three specific numerical examples (both two- and three dimensional) and study the effect of different noise measures.

used to model binary metal alloys [1], polymers [2] as well as cell proliferation and adhesion [3]. In material science, when a binary alloy is sufficiently cooled down, we observe a partial nucleation or spinodal decomposition, i.e., the material quickly becomes inhomogeneous. In fact, after a few seconds, material coarsening will be happened [4]. In polymer solutions and blends, the phase separation process is a dynamic process that one phase stable solution separates into two equilibrium phases upon changes in temperature, pressure, concentration, or even flow fields [5]. In these cases, the spinodal decomposition is described by the Cahn-Hilliard model [6].
The equation is a nonlinear partial differential equation of fourth-order in space and first order in time for which an analytical treatment is not possible. There are several numerical techniques to solve the equation including the finite element method (FEM) [7], isogeometric analysis based on finite element method [8], multigrid finite element [9], conservative nonlinear multigrid method [10], least squares spectral element method [11], Monte Carlo methods [12], radial basis functions (RBF) [13] and meshless local collocation methods [14]. A finite element error analysis of the equation is given in [15]. Adaptive finite elements can also be applied to solve the equation using residuals based a posteriori estimates [16,17].
A difficulty of the numerical analysis of the Cahn-Hilliard equation is the discretization of the fourth-order operator. Here, after converting the fourth-order equation into a system of two second-order equations (by introducing an auxiliary variable) and writing the variational formulation, the Ciarlet-Raviart mixed finite element method is used for the spatial discretization. The method has been implemented for the damped Boussinesq equation by the authors [18] and they considered the convergence rate and the stability for the semi-discretization scheme and the fully discretized method. For the Cahn-Hilliard equation, the technique was used in [19,20] for the space discretization.
The stochastic Cahn-Hilliard equation was first considered by Cook [21]. The system allows for considering thermal fluctuations directly in terms of the Cahn-Hilliard-Cook (CHC) equation by a conserved noise source term. The thermal fluctuations play an essential role in the early stage of phase dynamics such as initial dynamics of nucleation [22,23]. Some authors, such as Binder [24] and Pego [25], have expressed the belief that only the stochastic version can correctly describe the whole decomposition process in a binary alloy [26]. In [27], as another numerical approach, the authors employed the direct meshless local Petrov-Galerkin (DMLPG) to solve the stochastic Cahn-Hilliard-Cook and stochastic Swift-Hohenberg equations.
Multilevel Monte Carlo (MLMC) [28] is a simple and efficient computational technique to estimate the expected value of a random process. Using the method enables us to decrease the computational costs noticeably. The multilevel method was implemented to solve the stochastic elliptic equations, e.g., the drift-diffusion-Poisson system with uniformly distributed random variables [29] and quasi-random points [30]. In [31], the convergence and complexity of the MLMC using Galerkin discretizations in space and a Euler-Maruyama discretization in time for the parabolic equations were explained in details. The technique was used in [32] for solving parabolic (heat equation) and hyperbolic (advection equation) driven by additive Wiener noise Generally, for the time-dependent stochastic problems, the total error consists of the spatial error (due to the finite element method), the time discretization error (due to the Euler-Maruyama technique) and the statistical error (number of samples). We already know that for the space discretization, fine meshes are needed (specifically for the curved surfaces) which lead to the higher computational complexity. The multilevel Monte Carlo method uses hierarchies of meshes for time and space approximations in the sense that the number of samples and mesh sizes (as well as time steps) on the different levels are chosen such that the errors are equilibrated. For the stochastic Cahn-Hilliard-Cook equation, we strive to determine an optimal hierarchy of meshes, number of samples and time intervals which minimize the total computational work. As a result, we give a-priori estimates on the explained error contributions. In this paper, we use the MLMC-FEM for the fourth-order stochastic equations and calculate the mild solution of the Cahn-Hilliard-Cook equation. In fact, we estimate the total computational error according to the three error contributions. Then, we strive to minimize the computational complexity with respect to a given error tolerance. This procedure is compared to the Monte Carlo method.
The rest of the paper is organized as follows. In Sect. 2, we explain the Cahn-Hilliard and the Cahn-Hilliard-Cook equations with their boundary conditions. Then, we describe how the Ciarlet-Raviart mixed finite element can be used to convert the stochastic equation to a system of second-order equations. In Sect. 3, we demonstrate the implementation of the MLMC-FEM for the time-dependent stochastic equations. In Sect. 4, we give three numerical examples according to two different initial conditions. The solutions of the stochastic equation (the concentration) and the optimization (the optimal hierarchies) are given in this section. Finally, the conclusions are drawn in Sect. 5.
with the Neumann boundary conditions We consider the initial condition at t = 0 as where ν denotes the unit outward normal of the boundary and is a bounded domain in R d (d = 1, 2, 3). The solution u is a rescaled density of atoms or concentration of one of the material components where, in the most applications u ∈ [−1, 1]. We should note that M is the mobility (here a constant) and the variable is a positive constant. The equation arises from the Ginzburg-Landau free energy The above free energy includes the bulk energy F(u) and the interfacial energy (the second term). A popular example of a nonlinear function is The Cahn-Hilliard-Cook equation presents a more realistic model including the internal thermal fluctuations. It can be derived from (1) by adding the thermal noise as where ξ indicates the colored noise (here white noise) and σ is the noise intensify measure.

Ciarlet-Raviart mixed finite element
To construct a mixed finite element approximation of the Cahn-Hilliard-Cook equation, we first find its weak formulation. For this purpose, we define the auxiliary variable Therefore, the Cahn-Hilliard-Cook equation can be rewritten in the form The weak formulation of (8) is given by seeking where Now let τ h be a family of triangulations of into a finite number of elements (simplex) such that We assume that each element has at least one face on ∂ and k 1 , k 2 ∈ τ h have only one common vertex or a whole edge. Now we define and P n is the space of all polynomials of degree at most n ≥ 1. The semi-discrete Galerkin approximation of the solutions (9a)-(9b) may be defined as a pair of approximations hold.

Full discretization scheme
In this section we apply a fully discretize scheme based the mild solution of (8). In order to obtain the fully discretized scheme, we first rewrite the variational formulation of (9) as follows: The mixed finite element formulation of (15) is defined by Now we can rewrite (8) in the following abstract evolution equation where A is the negative Neumann Laplacian considered as an unbounded operator in the Hilbert space H = L 2 ( ), which is the generator of an analytic semigroup (S(t), t ≥ 0) on H [33]. The initial value X 0 is deterministic and W is a cylindrical Wiener process in H (i.e., the spatial derivative of a space-time white noise) with respect to a filtered probability Here, β j,k j,k∈N indicates a family of real-valued, identically distributed independent Brownian motions and μ j,k j,k∈N denote the eigenvalues (here, μ j,k = 1 since W (t) is cylindrical) [34]. Therefore, the Cahn-Hilliard-Cook equation has a continuous mild solution where t ∈ [0, T ], X : [0, T ] × → H and S(t) = e −t A 2 used as the analytic semigroup generated by −A 2 . The existence of the mild solution X was shown in [35]. Considering where C is a constant which depends on T . Also, for 0 ≤ s < t ≤ T , there exists a constant C(T ) such that the mild solution satisfies the inequality [31] In order to estimate the mild solution we use finite elements for space discretization and the semi-implicit Euler-Maruyama scheme in time direction. Let us assume that V ( ∈ N 0 ) is a nested family of finite element subsequences of H with refinement level > 0 and refinement size h ( ∈ N 0 ). Defining the analytic semigroup S = e −t A 2 , for t ∈ T , the semidiscrete problem (20) has the form For the time direction, we approximate the time discretization with step sizes δt ζ = T r −ζ where r > 1. Therefore, for ζ ∈ N 0 , we define the sequence of equidistant time discretization. In the computational geometry ( ), we estimate the mild solution X , with a finite element discretization. In other words, we suppose that the domain can be partitioned into quasi-uniform triangles or tetrahedra such that sequences {τ h } ∞ =0 of regular meshes are obtained. For any ≥ 0, we denote the mesh size of τ h by Uniform refinement of the mesh can be achieved by regular subdivision. This results in the mesh sizes where h 0 denotes the mesh size of the coarsest triangulation and r > 1 is independent of .

Multilevel Monte Carlo finite element method
The Monte Carlo method is a simple and efficient computational technique to solve SPDEs. As already mentioned, we use Euler-Maruyama to solve the equation on [0, T ] and the finite element method for the space discretization.
In order to obtain the mean square error (MSE) of ε, we require δt = O(ε) (for the time discretization). The Monte Carlo error (statistical error . Using a finite element scheme also gives rise to O(ε −d/α ), where α is the convergence rate of the discretization error. Therefore, by taking M samples, T /δt time steps and h as the mesh size, we have the following total cost It is obvious that for high dimensional geometries (i.e., d = 2, 3), the computational cost increases noticeably. Multilevel Monte Carlo finite element method (MLMC-FEM) is an efficient alternative to the Monte Carlo method to decrease the cost. In the time discretization, the general idea of the technique is using a hierarchy of the time steps, i.e., δt ( ∈ N 0 ) at different levels (instead of a fixed time step). For the space discretization, we use the mesh refinement (25) to obtain the mesh size at level which leads to In this section we strive to estimate the expectation of the mild solution on level L. First, for a given Hilbert space For can be defined as where taking expected value of the above equality leads to In order to approximate E[Y − Y −1 ] we can use the Monte Carlo estimator E M [Y − Y −1 ] (i.e., expectation of the difference of Y and Y −1 ) with independent number of samples M at level . Therefore, (31) can be estimated as In this part, we first provide the error bound for the single level Monte Carlo finite element. Then, using the obtained results, we achieve the error bound of the multilevel Monte Carlo considering the principal discretization error, i.e., spatial discretization (using finite element method), time stepping errors (due to the Euler-Maruyama technique) and statistical (sampling) error.
Lemma 1 [29] For any number of samples M ∈ N and for Y ∈ L 2 ( ; H ), the inequality According to Lemma 1 for , ζ ∈ N 0 and t ∈ ζ , we have the inequality where X ,ζ (t) is the discrete mild solution at level and time interval ζ . In order to estimate the discretization error which stems from the spatial discretization and time stepping we define the following lemma.

Lemma 2 [31]
Let X be the solution of (20) and X ,ζ be the sequence of discrete mild solution (i.e., the solution of (23)).

Then, there is a constant C(T ) such that for all
Hence, the total computational error is given by [31] sup In order to prove, we add and subtract the term E[X ,ζ (t)] to the left side and use the triangle inequality, Lemma 1, Lemma 2 and (34) to obtain the error bound. Now we couple the space and time errors and choose δt ζ h 2 (i.e., δt ζ = O(h 2 ) and O(δt ζ ) = h 2 ). Therefore, the total work W for the given spatial discretization level ∈ N 0 is estimated by After estimating the error bounds for single level Monte Carlo, we provide the multilevel Monte Carlo error bounds. By using ζ = 2 (due to δt ζ h 2 ), for ∈ N 0 we consider h r − and define δt := T r −2 . Therefore, for the full discretization in the multilevel setting, we redefine (24) as and we define the multilevel Monte Carlo estimator The computational error is given by For fixed L ∈ N and any chosen t L k(L) ∈ L , we split the error into two parts, i.e., discretization error and statistical error (see Theorem 4.5 in [31]), to obtain Now we make the following convergence assumptions for the prescribed errors. The below assumption is used to estimate the convergence rate of the discretization error Regarding the statistical error, the next assumptions are made. Hence, the total error can be estimated as The total work can be estimated by summing up the work of each level, i.e.,  Table 2 The optimal values of MC-FEM with respect to different prescribed errors Now we define an optimization problem which minimizes the computational work (46) such that the error is less or equal than a given error tolerance (ε). In other words, for estimating the mild solution on level L, we estimate the optimal hierarchies of (h , M , L) =L =0 such that In the problem we have M > 1, h 0 > 0 and r > 1. The exponent (α) as well as the coefficients (C 1 , C 2 , C 3 ) must be determined analogously. Finally, we should note that for Monte Carlo method the optimization problem with respect to (26) can be written as where again the optimization problem is over M > 1 and h > 0.

Numerical results
In

A 2D example
As the first example, we take u 0 = 0.25 + 0.1ω as the initial condition. The random variable ω is uniformly distributed between 0 and 1. The computational geometry ( ) is a circle with radius r = 2 and zero center point. As the first step, we try to solve the optimization problems (for MLMC-FEM and MC-FEM). This will enable us to find the optimal number of samples and mesh sizes. As explained in Sect. 2.1, we use the Ciarlet-Raviart mixed finite element with P1 elements. The estimation of the exponent α is crucial, however, it relates to the order of polynomials. Due to the fact that the exact solution of the stochastic equation is not available, we calculate the error between different mesh sizes and h = 0.01 (as the reference solution) at T = 5. Figure 1 depicts α with respect to different mesh sizes (here uniform refinement). The simulations show α = 0.952, C 1 = 0.51, where the exponent agrees very well with the order of P1 finite element technique (linear elements). The rest of the coefficients is estimated as C 2 = 0.066, C 3 = 0.223. Now we are ready to solve the optimization problem (47) with respect to the aforementioned parameters. In order to solve the optimization problem, we apply interior method where the details of the technique are given in [29]. The optimal hierarchies of the MLMC-FEM are shown in Table 1.
As the next step, in order to compare the efficiency of the multilevel setting with the Monte Carlo simulation, we solve the optimization given in (48). Again, the optimal mesh size and the optimal number of samples are given in Table 2 where the same convergence rate (α) is used. Finally, we draw a comparison between MLMC-FEM and MC-FEM which is shown in Fig. 2 Fig. 3. It is shown that from T = 0.2 to T = 1, a slow coarsening happens. Here, we use ε = 0.05 in the sense that the solution at the last level (here L = 2) is shown in the figure (see Table 1 for the optimal mesh size and number of samples). In order to study the noise effect we solve the deterministic equation with the same mesh size as illustrated in Fig. 4.

The 3D examples
Here we choose a more complicated example and use MLMC-FEM and CR-MFE to obtain the solution (expected value) of CHC equation in a cubic geometry, i.e., = where (x, y, z) is a point on the cube. The same procedure for solving the optimization problem can be used, however, due to the three-dimensional structure we set d = 3. We should note that as ζ = 2 , we define the optimal time interval as δt ζ h 2 . First, we consider the effect of the noise intensify measure in the sense that the deterministic case (σ = 0) is compared with the stochastic equation (σ = 0.1, 0.2, 0.3).
The results are shown in Fig. 5 at T = 1. Clearly, higher σ affects the concentration mostly. Similar to the 2D example, we consider the effect of time on the concentration. The results are shown in Fig. 6 for different times from T = 0.1 to T = 10. It illustrates that the initial homogeneous phase quickly segregates (at T = 0.1), however, after the segregation the domain starts to slowly coarsen in time.
In the next step, we use the Monte Carlo finite element method to compare the effect of the number of grids. Here, two mesh sizes, i.e., h = 0.5 (with 137 3 nodes) and h = 0.1 (with 6651 3 nodes) are employed and the results are shown in Figs. 7 and 8 for stochastic and deterministic cases, respectively. We solved the CHC equation with σ = 0.15 and compared its solution with the deterministic case (σ = 0) at time T = 5. It shows that the mesh size does not considerably affect the solution.
Finally, we consider the second 3D example (a torus). The first comparison is regarding the evolution of the concentration which is illustrated in Fig. 9 where in the simulations σ = 0.1 is used. In the second case, we study the effect of different noise measures from deterministic case to stochastic case with σ = 0.5 at T = 1. Here, the results are shown in Fig. 10. The simulations point out that the effect of σ = 0.4 and σ = 0.5 on the concentrations are more tangible.

Conclusions
In this paper, we considered the Cahn-Hilliard and Cahn-Hilliard-Cook equations as forth-order time-dependent equations. As the first step, after defining an auxiliary variable, we converted the equation into a system of second-order time-dependent problems. Then, we presented a variational formulation for the system and used the Ciarlet-Raviart We have already shown that for the stochastic timedependent problems, the computational cost of the Monte Carlo finite element method is O(ε −(2+1+d/α) ). Applying the multilevel technique for this problem reduces noticeably the computational costs. In a two-dimensional problem, the opti-mal hierarchies h , r , M , δt =L =0 reduce the complexity to O(ε −3.27 ) as certified in numerical example.
We showed three numerical examples with two specific initial conditions. We estimated the solution of stochastic/deterministic Cahn-Hilliard equation for different time intervals. As a result, we demonstrated distinctive coarsening and phase separation. For the stochastic equation, we studied the effect of noise measure, for showing that more σ intensifies the noisy concentration.