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.

Introduction

The Cahn–Hilliard equation is a robust mathematical model for describing different phase separation phenomena, from co-polymer systems to lipid membranes. The equation is 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.

Cahn–Hilliard–Cook equation

J. W. Cahn and J. E. Hilliard proposed the Cahn–Hilliard (CH) equation. The equation is a mathematical physics model that describes the process of phase separation. The CH equation is as followswith the Neumann boundary conditionsWe consider the initial condition at aswhere denotes the unit outward normal of the boundary and is a bounded domain in . The solution u is a rescaled density of atoms or concentration of one of the material components where, in the most applications . 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 energyThe above free energy includes the bulk energy F(u) and the interfacial energy (the second term). A popular example of a nonlinear function isThe 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 variableTherefore, the Cahn–Hilliard–Cook equation can be rewritten in the form The weak formulation of (8) is given by seeking such that whereNow let be a family of triangulations of into a finite number of elements (simplex) such thatWe assume that each element has at least one face on and have only one common vertex or a whole edge. Now we defineand is the space of all polynomials of degree at most . The semi-discrete Galerkin approximation of the solutions (9a)–(9b) may be defined as a pair of approximations for which the equalities 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:

Find such that The mixed finite element formulation of (15) is defined by such that Now we can rewrite (8) in the following abstract evolution equationwhere A is the negative Neumann Laplacian considered as an unbounded operator in the Hilbert space , which is the generator of an analytic semigroup on H [33]. The initial value 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 space defined asHere, indicates a family of real-valued, identically distributed independent Brownian motions and denote the eigenvalues (here, since W(t) is cylindrical) [34]. Therefore, the Cahn–Hilliard–Cook equation has a continuous mild solutionwhere , and used as the analytic semigroup generated by . The existence of the mild solution X was shown in [35]. Considering , for all the solution X satisfies [31]where C is a constant which depends on T. Also, for , 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 () is a nested family of finite element subsequences of H with refinement level and refinement size ). Defining the analytic semigroup , for , the semidiscrete problem (20) has the formFor the time direction, we approximate the time discretization with step sizes where . Therefore, for , we define the sequenceof 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 of regular meshes are obtained. For any , we denote the mesh size of byUniform refinement of the mesh can be achieved by regular subdivision. This results in the mesh sizeswhere denotes the mesh size of the coarsest triangulation and 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 (for the time discretization). The Monte Carlo error (statistical error) is (where M is the number of samples) which yields . Using a finite element scheme also gives rise to , where is the convergence rate of the discretization error. Therefore, by taking M samples, time steps and h as the mesh size, we have the following total costIt is obvious that for high dimensional geometries (i.e., ), 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., () 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 toIn this section we strive to estimate the expectation of the mild solution on level L. First, for a given Hilbert space the space, is defined to be the space of all measurable functions such thatFor the standard Monte Carlo estimator can be defined aswhere for , indicated a sequence of i.i.d. copies of . Let be a sequence of random variables such that , we can write astaking expected value of the above equality leads toIn order to approximate we can use the Monte Carlo estimator (i.e., expectation of the difference of and ) with independent number of samples at level . Therefore, (31) can be estimated asIn 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 and for , the inequalityholds for the MC error, where .

According to Lemma 1 for and , we have the inequalitywhere 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 be the sequence of discrete mild solution (i.e., the solution of (23)). Then, there is a constant C(T) such that for all , we have

Hence, the total computational error is given by [31]In order to prove, we add and subtract the term 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 (i.e., and ). Therefore, the total work for the given spatial discretization level is estimated byAfter estimating the error bounds for single level Monte Carlo, we provide the multilevel Monte Carlo error bounds. By using (due to ), for we consider and define . Therefore, for the full discretization in the multilevel setting, we redefine (24) asand we define the multilevel Monte Carlo estimatorThe computational error is given byFor fixed and any chosen , we split the error into two parts, i.e., discretization error and statistical error (see Theorem 4.5 in [31]), to obtainNow 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 assumptionsare made. Hence, the total error can be estimated asThe total work can be estimated by summing up the work of each level, i.e.,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 such thatIn the problem we have and . The exponent () as well as the coefficients () must be determined analogously. Finally, we should note that for Monte Carlo method the optimization problem with respect to (26) can be written aswhere again the optimization problem is over and .

Numerical results

In this section, we present three numerical examples of the stochastic Cahn–Hilliard–Cook equations where in all cases the optimal MLMC-FEM is used to obtain the solution. Due to the fact that the examples are real-world problems, their exact solutions are not given. The simulations are performed using MATLAB 2017b software on an Intel Core i7 machine with of memory. In all examples, is used and the constant mobility is applied. For the nonlinear term (i.e., ), we use Newton’s method where several iterations are needed to reach the stopping tolerance (here ). In each iteration, the built-in direct solver is employed to solve the linearized system.

A 2D example

As the first example, we take as the initial condition. The random variable is uniformly distributed between 0 and 1. The computational geometry () is a circle with radius 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 (as the reference solution) at . Figure 1 depicts with respect to different mesh sizes (here uniform refinement). The simulations show , , where the exponent agrees very well with the order of P1 finite element technique (linear elements). The rest of the coefficients is estimated as , . 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.

Fig. 1

The discretization error as a function of different mesh sizes where as the convergence rate is obtained

Table 1

The optimal hierarchies of MLMC-FEM with respect to different prescribed errors

ε	h_0	r	M_0	M_1	M_2	M_3	M_4
0.100	0.764	1.458	73	50	15	–	–
0.050	0.615	1.568	280	138	31	–	–
0.020	0.461	1.726	1672	524	85	–	–
0.010	0.370	1.856	6474	1445	184	–	–
0.005	0.580	1.990	187,700	46,448	4615	459	–

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. The results indicate that the multilevel method costs approximately and the computational work of Monte Carlo sampling is . The comparison indicates noticeably the efficiency of the MLMC-FEM.

Table 2

The optimal values of MC-FEM with respect to different prescribed errors

ε	0.1	0.05	0.02	0.01	0.005
h	0.108	0.052	0.020	0.010	0.005
M	3	12	73	289	1 152

Fig. 2

A comparison between the optimal work of MLMC-FEM and MC-FEM showing the efficiency of the multilevel technique is pronounced

Finally, we compare the evolution of the concentration at different times (from to ) (with ) where the obtained results are depicted in Fig. 3. It is shown that from to , a slow coarsening happens. Here, we use in the sense that the solution at the last level (here ) 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.

Fig. 3

The solution of the stochastic CHC equation at (top left), (top right), (bottom left) and (bottom right)

Fig. 4

The solution of the Cahn–Hilliard equation () at (top left), (top right), (bottom left) and (bottom right)

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., . The initial condition iswhere (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 . We should note that as , we define the optimal time interval as . First, we consider the effect of the noise intensify measure in the sense that the deterministic case () is compared with the stochastic equation (). The results are shown in Fig. 5 at . 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 to . It illustrates that the initial homogeneous phase quickly segregates (at ), however, after the segregation the domain starts to slowly coarsen in time.

Fig. 5

A comparison between different noise intensify measure, i.e., (top left), (top right), (bottom left) and (bottom right) at

Fig. 6

The expected value of the solution of CHC equation at (top left), (top right), (bottom left) and (bottom right). Here, and

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., (with nodes) and (with 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 and compared its solution with the deterministic case () at time . It shows that the mesh size does not considerably affect the solution.

Fig. 7

A comparison between two mesh sizes (left) and (right) at for the stochastic Cahn–Hilliard–Cook equation

Fig. 8

A comparison between two mesh sizes, i.e., (left) and (right) at for the deterministic Cahn–Hilliard equation

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 is used. In the second case, we study the effect of different noise measures from deterministic case to stochastic case with at . Here, the results are shown in Fig. 10. The simulations point out that the effect of and on the concentrations are more tangible.

Fig. 9

The evolution of the solution of the times (top left), (top right), (bottom left) and (bottom right)

Fig. 10

A comparison between different noise intensify measure, i.e., (top left), (top right), (middle left) and (middle right), (bottom left) and (bottom right) at

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 mixed finite element method. Afterwards, we rewrote the equation as a stochastic ODE in order to estimate its mild solution u(t).

We have already shown that for the stochastic time-dependent problems, the computational cost of the Monte Carlo finite element method is . Applying the multilevel technique for this problem reduces noticeably the computational costs. In a two-dimensional problem, the optimal hierarchies reduce the complexity to 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.