Remarks on numerical experiments of Allen-Cahnequations with constraint via Yosida approximation

. We consider a one-dimensional Allen{Cahn equation with constraint from the view-point of numerical analysis. Our constraint is the subdiﬀerential of the indicator function on the closed interval, which is the multivalued function. Therefore, it is very diﬃcult to make numerical experiments of our equation. In this paper we approximate our constraint by Yosida approximation. Then, we study the approximating system of our original model numerically. In particular, we give the criteria for the standard forward Euler method to give stable numerical experiments of our approximating equation. Moreover, we give some numerical experiments of approximating equation.


Introduction
In this paper, for each ε ∈ (0, 1] we consider the following Allen-Cahn equation with constraint from the view-point of numerical analysis: where 0 < T < +∞ and u ε 0 is a given initial data. Also, ∂I [−1,1] (·) is the subdifferential of the indicator function I [−1,1]  More precisely, ∂I [−1,1] (·) is a set-valued mapping defined by (1.5) The Allen-Cahn equation was proposed to describe the macroscopic motion of phase boundaries. In the physical context, the function u ε = u ε (t, x) in (P) ε :={(1.1), (1.2), (1.3)} is the nonconserved order parameter that characterizes the physical structure. For instance, let v = v(t, x) be the local ratio of the volume of pure liquid relative to that of pure solid at time t and position x ∈ (0, 1), defined by where B r (x) is the ball in R with center x and radius r and |B r (x)| denotes its volume.
Note that the constraint ∂I [−1,1] (·) is the multivalued function. Therefore, it is very difficult to make numerical experiments of (P) ε . Recently, Farshbaf-Shaker et. al [11] gave the results of the limit of a solution u ε and an element of ∂I [−1,1] (u ε ), called the Lagrange multiplier, to (P) ε as ε → 0. Moreover, Farshbaf-Shaker et. al [12] gave numerical experiments to (P) ε via the Lagrange multiplier in one dimension of space for sufficient small ε ∈ (0, 1]. Also, they considered some approximating method. In fact, for δ > 0, they use the following Yosida approximation (∂I [−1,1] ) δ (·) of ∂I [−1,1] (·) defined by: where [z] + is the positive part of z. For each δ > 0, they considered the following approximation problem of (P) ε : x ∈ (0, 1). Then, they gave the following numerical result to (P) ε δ by the standard explicit finite difference scheme to (P) ε δ (see [12,Remark 5.3]): From Figure 1, we easily see that we have to choose the suitable constants ε, δ and mesh size of time △t and space △x in order to get stable numerical results of (P) ε δ . So, in this paper, for each ε > 0 and δ > 0, we give the criteria for the standard explicit finite difference scheme to give stable numerical experiments of (P) ε δ . To this end, we first consider the following ODE problem, denoted by (E) ε δ : Then, we give the criteria to get stable numerical experiments of (E) ε δ . Also, we give some numerical experiments of (E) ε δ . Moreover, we show the criteria to get stable numerical experiments of PDE problem (P) ε δ . Therefore, the main novelties found in this paper are the following: (a) We give the criteria to give stable numerical experiments of the ODE problem (E) ε δ . Also, we give numerical experiments to (E) ε δ for sufficient small ε ∈ (0, 1].

(b)
We give the criteria to give stable numerical experiments of the PDE problem (P) ε δ . Also, we give numerical experiments to (P) ε δ for sufficient small ε ∈ (0, 1].
The plan of this paper is as follows. In Section 2, we recall the solvability and convergence result of (E) ε δ . In Section 3, we consider (E) ε δ numerically. Then, we prove the main result (Theorem 3.1) corresponding to the item (a) listed above. Also, we give numerical experiments to (E) ε δ for sufficient small ε ∈ (0, 1] and δ ∈ (0, 1]. In Section 4, we recall the solvability and convergence result of (P) ε δ . In the final Section 5, we consider (P) ε δ from the view-point of numerical analysis. Then, we prove the main result (Theorem 5.1) corresponding to the item (b) listed above. Also, we give numerical experiments to (P) ε δ for sufficient small ε ∈ (0, 1] and δ ∈ (0, 1].

Notations and basic assumptions
Throughout this paper, we put H := L 2 (0, 1) with usual real Hilbert space structure, and denote by (·, ·) H the inner product in H. Also, we put V := H 1 (0, 1) with the usual norm In Sections 2 and 4, we use some techniques of proper (that is, not identically equal to infinity), l.s.c. (lower semi-continuous), convex functions and their subdifferentials, which are useful in the systematic study of variational inequalities. So, let us outline some notations and definitions. Let W be the real Hilbert space with the inner product (·, ·) W . For a proper, l.s.c. and convex function ψ : The subdifferential of ψ is a possibly multi-valued operator in W and is defined by For various properties and related notions of the proper, l.s.c., convex function ψ and its subdifferential ∂ψ, we refer to a monograph by Brézis [4].
, if the following conditions are satisfied: (ii) The following equation holds: We easily see that the problem (E) ε δ can be rewritten as in an abstract framework of the form: Therefore, applying the Lipschitz perturbation theory of abstract evolution equations (cf. [5,14,21]), we can show the existence of a solution u ε δ to (EP) ε δ on [0, T ] for each ε ∈ (0, 1] and δ ∈ (0, 1] in the sense of Definition 2.1. Thus, the proof of Proposition 2.1 has been completed. Next, we recall the convergence result of (E) ε δ as δ → 0. To this end, we recall a notion of convergence for convex functions, developed by Mosco [17]. Definition 2.2 (cf. [17]). Let ψ, ψ n (n ∈ N) be proper, l.s.c. and convex functions on a Hilbert space W . Then, we say that ψ n converges to ψ on W in the sense of Mosco [17] as n → ∞, if the following two conditions are satisfied: It is well known that the following lemma holds. Therefore, we omit the detailed proof.  [17] (2.3) as δ → 0.
By Lemma 2.1 and the general convergence theory of evolution equations, we easily get the following result.
and u ε is the unique solution of the following problem (E) ε on [0, T ]:

Stable criteria and numerical experiments for (E) ε δ
In this Section we consider (E) ε δ from the view-point of numerical analysis. Remark 3.1. Note from Proposition 2.2 that (E) ε δ is the approximating problem of (E) ε . Also note from (1.5) that the constraint ∂I [−1,1] (·) is the multivalued function. Therefore, it is very difficult to study (E) ε mumerically.
In order to make numerical experiments of (E) ε δ via the standard forward Eular method, we consider the following explicit finite difference scheme to (E) ε δ , denoted by of (DE) ε δ : where △t is the mesh size of time and N t is the integer part of number T /△t. We easily see that u n is the approximating solution of (E) ε δ at the time t = n△t Clearly, the explicit finite difference scheme (DE) ε δ converges to (E) ε δ as △t → 0 since (DE) ε δ is the standard time discretization scheme for (E) ε δ . Here, we give the unstable numerical experiment of (DE) ε δ in the case when T = 0.002, ε = 0.003, δ = 0.01, the initial data u ε 0 = 0.1 and the mesh size of time △t = 0.000001: From Figure 2, we easily see that we have to choose the suitable constants ε, δ and mesh size of time △t in order to get stable numerical results of (DE) ε δ . Now, let us mention our first main result in this paper, which is concerned with the criteria to give stable numerical experiments of (DE) ε δ .
Also, by (3.2) and (3.5) we observe that f δ (u n ) ≤ 0 for all n ≥ 0. Therefore, we observe from (3.3) that Therefore, we infer from (3.5) and (3.8) that {u n ; n ≥ 0} is a bounded and increasing sequence with respect to n. Thus, there exist a subsequence {n k } of {n} and a point u ∞ ∈ R such that n k → ∞ as k → ∞ and By taking the limit in (3.3), we easily observe from the continuity of f δ (·) that u ∞ = 1/(1−δ), which is the zero point of f δ (·). Hence, taking into account of the uniqueness of the limit point, the proof of (i) has been completed. Next, we show (ii). To this end, we assume that △t ∈ ( ) . Then, we can find the minimal number n 0 ∈ N so that Taking into account of (3.11), u 0 ∈ (0, 1] and 1 + △t we can find the minimal number n 0 ∈ N so that u n0 > 1 and u i ∈ (0, 1] for all i = 0, 1, · · · , n 0 − 1. Also, by (3.4) we observe that thus, (3.10) holds.
To show (ii), we put Then, we observe from (3.2) and (3.3) that Therefore, we observe from (3.13) and τ ∈ (1, 2) that the zero point is in the interval between u n0 and u n0+1 . Also, by (3.12) we observe that which implies that Therefore, by (3.10), (3.13) and by repeating the procedure as above, we observe that for all n ≥ n 0 (3.14) and u n oscillates around the zero point 1/(1 − δ) for all n ≥ n 0 . Also, we observe from (3.12) and (3.14) that Therefore, by τ ∈ (1, 2), (3.14) and (3.15), there exist a subsequence {n k } of {n} such that u n k oscillates and converges to 1/(1 − δ) as k → ∞. Hence, taking into account of the uniqueness of the limit point, the proof of (ii) has been completed.

Remark 3.3. By (3.3) we easily see that
By (ii) of Theorem 3.1, we observe that u n oscillates and converges to the zero point of f δ (·) in the case when △t ∈ ( . However, in the case when △t = 2δε 2 /(1 − δ), we have the following special case that the solution to (DE) ε δ does not oscillate and coincides with the zero point of f δ (·) after some finite number of iteration.
Then, we easily observe that:

The case when △t = 0.000001
Now we consider the case when △t = 0.000001. In this case, we have: which implies that (i) of Theorem 3.1 holds. Thus, we have the following stable numerical result of (DE) ε δ . Namely, the solution to (DE) ε δ converges to the stationary

The case when △t = 0.000002
Now we consider the case when △t = 0.000002. In this case, we have: which implies that (ii) of Theorem 3.1 holds. Thus, we have the following numerical result of (DE) ε δ that the solution to (DE) ε δ oscillates and converges to the stationary  Figure 4:

The case when △t = 0.000005
Now we consider the case when △t = 0.000005. In this case, we have: Therefore, we observe Remark 3.2. In fact, we have the following numerical result of (DE) ε δ that the solution to (DE) ε δ oscillates.   (1 − δ). In this case, we observe Remark 3.2. In fact, we have the following numerical result of (DE) ε δ that the solution to (DE) ε δ oscillates between three zero points of f δ (·).

Numerical result of Corollary 3.1
In this subsection, we consider Corollary 3.1 numerically. To this end, we use the follwing initial data: Then, we have the following numerical experiment of (DE) ε δ that Corollary 3.1 holds. Namely, we observe that (3.16) holds with n = 6:  Figure 8:

Conclusion of ODE problem (DE) ε δ
By Theorem 3.1 and numerical experiments as above, we conclude that (i) the mesh size of time △t must be smaller than δε 2 /(1 − δ) in order to get the stable numerical experiments of (DE) ε δ . (ii) we have the stable numerical experiments of (DE) ε δ with the initial data u ε 0 := (1 − δ) n−1 /(1 + δ) n , even if the mesh size of time △t is equal to 2δε 2 /(1 − δ).
The following variational identity holds: for all z ∈ V and a.e. t ∈ (0, T ).
Now, let us recall the solvability result of (P) ε δ on [0, T ]. Proposition 4.1. Let ε ∈ (0, 1] and δ ∈ (0, 1]. Assume the following condition: Also, we can show the existence of solutions to (P) ε δ on [0, T ] applying the abstract theory of evolution equations governed by subdifferentials. In fact, we define a functional φ ε δ on H by is the function defined in (2.1). Clearly, φ ε δ is proper, l.s.c. and convex on H with the effective domain D(φ) = V .
We easily see that the problem (P) ε δ can be rewritten as in an abstract framework of the form: Therefore, applying the Lipschitz perturbation theory of abstract evolution equations (cf. [5,14,21]), we can show the existence of a solution u ε δ to (PP) ε δ , hence, (P) ε δ , on [0, T ] for each ε ∈ (0, 1] and δ ∈ (0, 1] in the sense of Definition 4.1. Thus, the proof of Proposition 4.1 has been completed. Next, we recall the convergence result of (P) ε δ as δ → 0. Taking account of Lemma 2.1 (cf. (2.3)), we easily observe that the following lemma holds.
Then, φ ε δ (·) −→ φ ε (·) on H in the sense of Mosco [17] as δ → 0. By Lemma 4.1 and the general convergence theory of evolution equations, we easily get the following result. Proof. We easily observe that the problem (P) ε can be rewritten as in an abstract framework of the form: Therefore, by Lemma 4.1 and the abstract convergence theory of evolution equations (cf. [2,15]), we observe that the solution u ε δ to (PP) ε δ converges to the unique solution u ε to (PP) ε on [0, T ] as δ → 0 in the sense of (4.2). Note that u ε (resp. u ε δ ) is the unique solution to (P) ε (resp. (P) ε δ ) on [0, T ] (cf. Proposition 4.1). Thus, we conclude that Proposition 4.2 holds.

Stable criteria and numerical experiments for (P) ε δ
In this Section we consider (P) ε δ from the view-point of numerical analysis. Remark 5.1. Note from Proposition 4.2 that (P) ε δ is the approximating problem of (P) ε . Also note from (1.5) that the constraint ∂I [−1,1] (·) is the multivalued function. Therefore, it is very difficult to study (P) ε numerically.
In order to make numerical experiments of (P) ε δ , we consider the following explicit finite difference scheme to (P) ε δ , denoted by of (DP) ε δ : ε 2 for n = 0, 1, 2, · · · , N t and k = 1, 2, · · · , N x − 1, u n 0 = u n 1 , u n Nx = u n Nx−1 for n = 1, 2, · · · , N t , where △t is the mesh size of time, △x is the mesh size of space, N t is the integer part of number T /△t, N x is the integer part of number 1/△x and x k := k△x. We easily see that u n k is the approximating solution of (P) ε δ at the time t n := n△t and the position x k := k△x.
From Figure 1, we easily see that we have to choose the suitable constants ε, δ, the mesh size of time △t and the mesh size of space △x in order to get stable numerical results of (DP) ε δ . Now, let us mention our second main result in this paper, which is concerned with the stability of (DE) ε δ .
By the similar arguments as above, we observe that the function is non-negative and continuous. In fact, it follows from (3.2) that the function attains a minimum value at z = −1. Therefore, we observe from (3.2) and (5.1) that . (5.11) Also, for any z ∈ [−1/(1 − δ), −1], we observe from (3.2) that Here we note from (5.1) that is non-decreasing and attains a minimum value at z = −1/ (1 − δ). Hence, we have: ] .

(5.13)
Thus, we observe from (5.11) and (5.13) that which implies from (5.3) and (5.10) that Taking into account of Neumann boundary condition, namely, by u n+1 and u n+1 Nx = u n+1 Nx−1 , we observe from (5.9) and (5.14) that which implies that (5.3) holds for i = n + 1. Therefore, we conclude from the mathematical induction that (5.2) holds. Hence, the proof of (i) of Theorem 5.1 has been completed. Next, we show (ii) by the standard arguments. Namely, we reformulate (DP) ε δ as in the following form: Here by taking into account of Neumann boundary condition and initial condition, namely, by u n 0 = u n 1 and u n Nx = u n Nx−1 for all n ≥ 0, we observe that (5.15) is reformulated as in the following form: where we put Noting from (3.2) that By using the matrix as above, we observe that (5.15) can be rewritten as in the following form: where we put By the general theory, we observe that the eigenvalue λ j of matrix A is given by: and (5.24) Therefore, we observe from (5.23)-(5.24) that max 1≤k≤Nx−1 |u n k | is increasing with respect to n in the case when max 1≤k≤Nx−1 |u n k | ≤ 1. However, if u n k / ∈ [−1, 1] for some k = 1, 2, · · · , N x − 1, it follows from (5.1) and (5.18) that Therefore, the sum of all components in k-th row of A + B is the following: 1] for some k = 2, 3, · · · , N x − 2. (5.26) Althought max 1≤k≤Nx−1 |u n k | is increasing with respect to n in the case when max 1≤k≤Nx−1 |u n k | ≤ 1 (cf. (5.23)-(5.24)), we conculde from (5.22) and (5.25)-(5.26) that (ii) of Theorem 5.1 holds. Thus, the proof of Theorem 5.1 has been completed.
Taking into account of Theorem 5.1, we give numerical experiments of (DP) ε δ as follows. To this end, we use the following numerical data: Numerical data of (DP) ε δ .
Also, we consider the following initial data u ε 0 (x) defined by Therefore, we have c 0 δε 2 1 − δ = 0.00001515151515 · · · > △t, which implies that the criteria condition (5.1) holds. Thus, we have the following stable numerical experiment of (DP) ε δ : Therefore, we have c 0 δε 2 1 − δ = 0.00000029696969 · · · < △t, which implies that the criteria condition (5.1) does not hold. Therefore, we have the following unstable numerical experiment of (DP) ε δ : which implies that the criteria condition (5.1) holds. Therefore, we have the following stable numerical experiment of (DP) ε δ : Remark 5.4. We observe from Theorem 5.1 that in order to get stable numerical results of (DP) ε δ , we have to choose the suitable constants ε, δ and the mesh size of time △t and space △x. Therefore, if we make a numerical experiment of (P) ε for sufficient small ε, we had better consider the original problem (P) ε by using a primal-dual active set method in [3], a Lagrange multiplier method in [12] and so on.

Conclusion of PDE problem (DP) ε δ
By Theorem 5.1 and numerical experiments as above, we conclude that the mesh size of time △t and space △x must be satisfied for some constant c 0 ∈ (0, 1), in order to get the stable numerical experiments of (DP) ε δ . Also, by Theorems 3.1 and 5.1, we conclude that the value δε 2 /(1 − δ) is very important to make numerical experiments of (DE) ε δ and (DP) ε δ .