STRUCTURE-PRESERVING FINITE DIFFERENCE SCHEMES FOR THE CAHN–HILLIARD EQUATION WITH DYNAMIC BOUNDARY CONDITIONS IN THE ONE-DIMENSIONAL CASE

. The structure-preserving ﬁnite diﬀerence schemes for the one di-mensional Cahn–Hilliard equation with dynamic boundary conditions are stu- died. A dynamic boundary condition is a sort of transmission condition that includes the time derivative, namely, it is itself a time evolution equation. The Cahn–Hilliard equation with dynamic boundary conditions is well-treated from various viewpoints. The standard type consists of a dynamic boundary condition for the order parameter, and the Neumann boundary condition for the chemical potential. Recently, Goldstein–Miranville–Schimperna proposed a new type of dynamic boundary condition for the Cahn–Hilliard equation. In this article, numerical schemes for the problem with these two kinds of dynamic boundary conditions are introduced. In addition, a mathematical result on the existence of a solution for the scheme with an error estimate is also obtained for the former boundary condition.


1.
Introduction. It is natural that volume conservation for the order parameter is required for the Cahn-Hilliard equation, because this system is well-known as the partial differential equation that describes the anti-diffusive phenomena of the solid-solid phase separation for a binary alloy. Therefore, the volume of the order parameter should be invariant with respect to time. If we use the Neumann boundary condition for the chemical potential, then the volume of the order parameter is conserved. The same is true under the standard type of dynamic boundary condition, more precisely, under a dynamic boundary condition for the order parameter and the Neumann boundary condition for the chemical potential given below as (DBC1). Recently, Goldstein-Miranville-Schimperna [17] proposed a new type of dynamic boundary condition ((DBC2) below). In this condition, the time derivative of the order parameter is combined with the normal derivative of the chemical potential. Roughly speaking, this problem represents a system between the Cahn-Hilliard equation in the bulk and on the boundary. In this case, the system has characteristic conservation structure called total mass conservation, that is, the sum of the volumes of the order parameter in the bulk and on the boundary is conserved. Therefore, from the viewpoint of numerical analysis, it is interesting to construct structure-preserving schemes related to each of these volume conservations.
Let T > 0 be a finite time and let L > 0 be the length of the one-dimensional material. In this article, we study the following one-dimensional Cahn-Hilliard equation: under two kinds of dynamic boundary condition. The unknowns u := u(x, t) and p := p(x, t) are the order parameter and the chemical potential, respectively, F is some derivative of the potential F , and γ is a positive constant. The first problem, under the standard dynamic boundary condition, is of the following form: which has been studied in many articles (see, e.g., [5,8,9,10,15,16,19]). The second problem, recently proposed by Goldstein-Miranville-Schimperna [17], is      ∂ t u(0, t) = ∂ x p(0, t), ∂ t u(L, t) = −∂ x p(L, t), p(0, t) = −γ∂ x u(0, t) + F (u(0, t)), p(L, t) = γ∂ x u(L, t) + F (u(L, t)), t ∈ (0, T ]. (DBC2) The Cahn-Hilliard equation (CH) with the condition (DBC2) is interpreted as the problem with non-permeable walls. In other words, this problem combines the Cahn-Hilliard equation (CH) in the bulk (0, L) and the Cahn-Hilliard equation (DBC2) on the boundary at x = 0 and x = L. We remark that the original problem was introduced in the three-dimensional case, where (DBC2) include the Laplace-Beltrami operator on the boundary, which plays the role of diffusion on the boundary. Related studies are found in [2] and [6] in the wider setting of a double-well potential F . In this article, we propose structure-preserving finite difference schemes for both (CH) with (DBC1) and (CH) with (DBC2). Suppose that the potential function F is in C 3 and is bounded from below. A well-known classical example of F is the polynomial double-well potential for constants q, r > 0. Let us define

1917
Then it is easily seen that the solution for (CH) with (DBC1) satisfies for all t ∈ (0, T ]. On the other hand, the solution for (CH) with (DBC2) satisfies for all t ∈ (0, T ]. Equations (1) and (3) represent total mass conservation, and equations (2) and (4) represent energy dissipation. In this article, structure-preserving means that these laws are inherited in the finite difference schemes in the discrete sense. We remark that the difference between these structure-preservations may result in large differences for the long-time behavior of each solution. Interesting related results can be found in [21] for (CH) with (DBC1) and [20] for a problem similar to (CH) with (DBC2). It is known that structure-preserving numerical schemes naturally assure the stability of numerical solutions. Moreover, from the perspective of mathematical analysis, various strategies for the continuous case such as the energy method can be applied to the discrete numerical scheme in a similar way (see, e.g., [22]). For a classical boundary condition such as the Neumann boundary condition or a Dirichlet boundary condition, systematic methods for deriving structure-preserving finite difference schemes are known. One example is the discrete variational derivative method (DVDM) proposed by Furihata [13]. However, these methods do not seem to apply directly to our problem because it is difficult to define a discrete version of the variational derivative for a conserved quantity with dynamic boundary conditions. We thus derive the schemes by trial and error. Here we remark that there are few numerical studies of (CH) with (DBC2) (see, e.g., [3,4,18] for finite element methods). In this sense, our numerical method is a finite-difference method, although it is challenging from the viewpoint of structure-preservation. Because volume conservation is an inherent in the setting of function spaces in these previous studies, this difficulty does not appear.
The rest of this article is organized as follows. In Section 2, we propose two numerical schemes. One is a scheme for (CH) with (DBC1), whose solution satisfies the discrete version of (1) and (2). The other is a scheme for (CH) with (DBC2), whose solution satisfies the discrete version of (3) and (4). We next compare the numerical solutions of the two schemes, in particular, we observe the behavior of the solutions by changing the parameter γ. The mathematical results in this article are gathered in Section 3. We obtain error estimates between the exact solution and the numerical solutions and discuss the existence and uniqueness of solutions for the numerical scheme in the case of (CH) with (DBC1).
2. Numerical schemes and simulations. In this section, we first present some notation and terminology. We then introduce our schemes, which inherit the discrete version of (1)- (4). Finally, we use our schemes to compare the numerical solutions for (CH) with (DBC1) and (DBC2).
The symbols ∂ x and ∂ t denote partial differential operators with respect to the space variable x and the time variable t, respectively. We also define the differential operators with respect to ξ, η, and γ by ∂ ξ , ∂ η , and ∂ γ respectively. In particular, in the case of a single variable function, we may also denote the derivatives by F , F and F . Let us define C m (Ω) as the function space of m-times continuously differentiable functions on Ω⊂ R. We remark that the domain Ω will be used in various situations, for instance, in some instances as a subset of [0, L], in others as the bounded ball {ξ | |ξ| ≤ R}, and so on. We also use the notation C m (Ω × [0, T ]), which is the function space of m-times continuously differentiable functions with respect to both the space and time variables. Let K and N be natural numbers, and define mesh sizes by ∆x := L/K and ∆t := T /N . In this article, we seek values at ((k − 1/2)∆x, n∆t) for k = 1, . . . , K and n = 0, 1, . . . , N . However, we formally define the numerical solution at not only k = 1, 2, . . . , K but also k = 0, K + 1. Observe that the points k = 0, K + 1 correspond to x = −∆x/2, L + ∆x/2 / ∈ [0, L]. We use the notation f (n) k for the value at ((k − 1/2)∆x, n∆t). We use bold font to denote (K + 1)-dimensional vectors with respect to the space variable, for example, k=0 and in particular f := (f k ) K k=0 in the single variable case. We remark that the bold font is used only for abbreviation of (K + 1)-dimensional vector although (K + 2)-dimensional vector will also appear. For approximations to derivatives and integrals, we follow the notation of [14]. That is, the difference Analogous to the expressions such as ∂ x f (a) = ∂ x f (x)| x=a , we denote the value of e.g. δ + n f (n) k at some point (m∆x, ∆t) by δ + n f ( ) m . For example, we will use and so on. In [14], which introduced DVDM, the trapezoidal rule was generally used as an approximation of the integral L 0 f (x)dx, and for various formulas and estimates for approximations such as summation by parts and the Sobolev inequality. In this article, we use the rectangle rule: Because the error for (7) is worse than for (6), the order of the error estimate for our problem may be also worse. On the other hand, it helps with extending the results to multi-dimensional problems. Indeed, in [14], the rectangle rule is used when the extension of DVDM to the multi-dimensional case is introduced. They only mention how to derive multi-dimensional structure-preserving schemes and give no mathematical treatment there. Thus we have to reconstruct the corresponding formulas for the rectangle rule, which is somewhat easier than for (6).

Figure 1. Relationship between x and k
To use our schemes for more general nonlinear cases, we define the difference quotient. Let Ω be a domain in R. For a function F ∈ C 1 (Ω) and ξ, η ∈ Ω, the difference quotient ∂F/∂(ξ, η) of F at (ξ, η) is defined as We denote the vector (∂F/∂(f k , g k )) K k=0 by ∂F/∂(f , g), and in the same fashion we denote (F (f k )) K k=0 by F (f ). Our scheme corresponding to (CH) with (DBC1) is given by for k = 1, 2, . . . , K, and the boundary condition corresponding to (DBC1) is

TAKESHI FUKAO, SHUJI YOSHIKAWA AND SAORI WADA
The boundary condition corresponding to (DBC2) is We define the discrete version of the volume M d in the bulk and the energy E d by The solution for (dCH) with (dDBC1) ((dCH) with (dDBC2), resp.) satisfies the following discrete version of the conservation laws corresponding to (1) and (2) ((3) and (4), resp.).
For the proof, we use the following lemma. This can easily be shown by direct calculation, and we thus omit the proof.
The following fundamental formula holds: (2) The following summation by parts formula holds: Substituting Proof of Theorem 2.1. We only show (10) because the other assertions can be proved in a similar manner. It follows from the summation by parts (14) that which completes the proof of (10).

Numerical computations.
This subsection is devoted to comparing the behavior of solutions for (CH) with (DBC1) and (DBC2) from our schemes. We define (9) and (10), the solution for (dCH) with (dDBC1) satisfies Similarly, if we let then the solution for (dCH) with (dDBC2) satisfies from (11) and (12). We choose K = 20 and N = 3000, and fix L = 1 so that ∆x = 1/20. However, we choose T depending on the situation. We fix the parameters q = 1 and r = 1, and thus consider the nonlinear function F (u) = u 4 /4 − u 2 /2. We also choose the initial value u(x, 0) = 0.1 sin(2πx) + 0.01 cos(4πx) + 0.06 sin(4πx) + 0.02 cos(10πx), which is the same as that used in [13]. All our numerical computations are simulated using Scilab. To obtain the next time-step for the solutions using our nonlinear scheme, we use the function "fsolve" in Scilab. We first check that the above quantities A (n) 2,d and M  From these numerical results, we make some conjectures related to the long-time behavior of the solutions. For (CH) with (DBC1), the target state may be the solution of some stationary problem with the following boundary conditions: That is, u is the solution of some Neumann boundary-valued problem. The target state for the boundary conditions should be for (CH) with (DBC2). This implies that the boundary-value of u solves some equation in F , that is, it depends on the equilibrium. From these simple pictures in the two cases, we find a difference in the long-time behavior of each of the solutions. We remark that asymptotic analysis as γ → 0 is also an interesting problem (see, e.g., [7,11,12]).

Existence of solutions and error estimates.
Here we give proofs of the existence of solutions for the schemes and error estimates. We restrict our argument For the problem with (dDBC2), we have not obtained any results yet. We will discuss the difficulty at the end of this section.
3.1. Preliminaries. Let us first define the discrete L ∞ d -norm and the discrete Dirichlet semi-norm D · for (K + 1)-dimensional vector by  Proof. Because Lf m = K k=1 f m ∆x, we see that It is easily seen that Thus we have Combining (16) with (17), we obtain Df . Therefore, In [1], the original Poincaré-Wirtinger inequality is introduced as We remark that, from the above proof, we can easily check that the corresponding discrete version |δ − k f k |∆x, m = 0, 1, . . . , K also holds. Next we define F (ξ, ξ; η, η) by which appears naturally by subtracting the difference quotients.

TAKESHI FUKAO, SHUJI YOSHIKAWA AND SAORI WADA
We also give several important properties of F . For proofs, we refer the reader to the second author's previous results in [22] and [23]. ( , then the following estimate holds: 3.2. Existence of solution. By applying the energy method proposed in [22], we have the following existence result.

Theorem 3.4 (Existence and uniqueness of solution). Let
for k = 1, 2, . . . , K, Once we show the existence of a unique fixed point for Ψ, we may regard the fixed point as {U can be explicitly written as by (23), it is sufficient to show the existence of a (K + 1)-dimensional vector U := where U k and U k satisfy (21)-(24), and V : We prepare an auxiliary mapping Φ : V → V . Given U (n) , we show that the mapping Φ is contraction into for R (n) := DU (n) . Obviously, showing the existence of unique fixed points for Ψ and Φ is equivalent. Observe that DU = DV and D U = D V . Moreover, it follows from (15), (21) and (24) that Assume that V ∈ X. From the definition of X, we see that M d (U ) = M d (U (n) ). By direct calculation, we have from (21)-(24)

1928
TAKESHI FUKAO, SHUJI YOSHIKAWA AND SAORI WADA because of (14) and the fact that ab ≤ a 2 + b 2 /4. Consequently, using the triangle inequality √ a 2 + b 2 ≤ |a| + |b|, we have From Proposition 3.1, we see that Next we define R (n) by Then, by Lemma 3.2, we see that U (n) , and by the fact that V ∈ X, we also have It follows from (19) that max k=1,2,...,K Thus we conclude that Next, we show that Φ is a contraction. Assume that for i = 1, 2 and k = 0, 1, . . . , K +1, and define U i : k=0 and { U i,k } K+1 k=0 (i = 1, 2) are connected by the relations for k = 1, 2, . . . , K, Subtracting these yields Thus we have

TAKESHI FUKAO, SHUJI YOSHIKAWA AND SAORI WADA
Because from Proposition 3.1 we have it follows from Lemma 3.1 that We observe that, by Lemma 3.3, Because U 1 − U 2 ∈ X, we see that M d (U 1 − U 2 ) = 0, and hence we obtain from Lemma 3.2 that and thus we conclude that then Φ is contraction into X, which completes the proof of local existence.
Step 2. From the energy conservation laws (Theorem 2.1), we have (9) and (10), we see that Thus we have and it follows from Lemma 3.2 that U (n) L ∞ d ≤ B 0 . Therefore, replacing R (n) and R (n) with B 0 and B 0 , respectively, we obtain the desired result. For more precise argument, we refer the reader to [22] and [23].
3.3. Error estimate. We next give an estimate of the error between the approximate solution and the exact solution. Observe that the approximate solution U given by  Proof. It follows from the definition of p and the boundary condition (DBC1) that From the same argument, we have Let us denote the errors e Subtracting the scheme from the original equation at t = (n + 1/2)∆t for k = 1, 2, . . . , K, we obtain where ζ 1,k , ζ 2,k are residue terms. It is easily seen from the same argument as [22] that ζ 1,k , ζ 2,k = O(∆t 2 ) + O(∆x 2 ) for k = 2, 3, . . . , K − 1 when ∆t, ∆x → 0. On the other hand, the convergence rates for k = 1, K are O(∆t 2 ) + O(∆x). Indeed, for e.g. k = 1 From the definition of u and the boundary condition we have Then we obtain Thus we conclude that ζ 2,1 = O(∆t 2 ) + O(∆x), and in the same fashion ζ 2,K = O(∆t 2 ) + O(∆x). Similarly, for ζ 1 we see that ζ 1,k = O(∆t 2 ) + O(∆x) at k = 1, K, since we can easily check that for k = 1, K.
On the boundary, thanks to (29) and (30), we have p,0 , For the boundary condition for u, because In the same fashion, we also obtain where ζ 3,0 , ζ 3,K = O(∆t 2 ) + O(∆x). Consequently, we arrive at the following expression for the error: p,K = 0. It follows from (14) that For the first term on the right-hand side, we use the estimate It follows from the Young inequality ab ≤ ( /2)a 2 + (1/2 )b 2 that By a calculation similar to that in the proof of existence using Lemmas 3.1 and 3.3 and Proposition 3.1, we obtain To estimate e  Then, because (a + b + c) 2 ≤ 2(a 2 + b 2 ) + ( /2)(a 2 + b 2 ) + (1 + 4/ ) c 2 , we conclude that From the same argument as the proof in [22], if we assume that ∆t < γ/A then we conclude that De (n) u ≤ C∆t 2 + C∆x, n = 0, 1, . . . , N.
Using the Poincaré-Wirtinger inequality again completes the proof.
3.4. Further remarks. We only discuss the mathematical results in the case (dDBC1) in spite of that we also introduce the numerical scheme in the case (dDBC2). The difficulty to treat the case (dDBC2) arises from several reasons. Comparing conservation laws (9) and (10) with (11) and (12), we see that boundary values which are not assured to be non-negative appear in the latter. It causes the difficulty to apply the existence proof of Theorem 3.4 to the problem with (dDBC2) directly. Moreover, according to the proof of Theorem 3.5, the boundary equations for error (31) play an important role. On the other hand, in the case (dDBC2) the residue terms ζ appear also in the boundary equation if we use the same argument as Theorem 3.5. Hence, now we have no idea to show the mathematical results for the problem with (dDBC2).
We also remark multi-dimensional problems. Our scheme are easily extended to multi-dimensional case formally by using the Voronoi mesh (see [14]). However, as we mentioned before, the term including Laplace-Beltrami operator appears on the boundary conditions in the multi-dimensional case even in (DBC1). Although it gives smoothing effect in the continuous case (see e.g. [5]), now we do not have idea how to extend our results to the problem with the Laplace-Beltrami operator.