# Efficient Implementation and Numerical Analysis of Finite Element Method for Fractional Allen-Cahn Equation

Authors: Guozhen Wang; Huanzhen Chen
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

We embed the fractional Allen-Cahn equation into a Galerkin variational framework and thus develop its corresponding finite element procedure and then prove rigorously its mathematical and physical properties for the finite element solution. Combining the merits of the conjugate gradient (CG) algorithm and the Toeplitz structure of the coefficient matrix, we design a fast CG for the linearized finite element scheme to reduce the computation cost and the storage toO(Mlog⁡M) and O(M), respectively. Numerical experiments confirm that the proposed fast CG algorithm recognizes accurately the mass and energy dissipation, the phase separation through a very clear coarse graining process, and the influences of different indices r of fractional Laplacian and different coefficients K,η on the width of the interfaces.


---

## Body

## 1. Introduction

As a typical phase-field model, the classical Allen-Cahn equation [1](1)∂tu-KΔu+ψ′u=0,was originally derived through the minimization of the Ginzburg-Landau free energy functional(2)∫Ω12∇u2+ψudx to describe the motion of antiphase boundaries in crystalline solids with the double-well potential ψ(u)=1/4η2(u2-1)2. Since then, the Allen-Cahn equation has been widely applied to many complicated moving interface problems, for example, vesicle membranes, nucleation of solids, and mixture of two incompressible fluids, etc. (cf. [2–8]), and many and many research results on its theories, applications, and numerics of the Allen-Cahn equation have been achieved; see the reviews [9–15] and the references cited therein.To recognize the influences of the long-range interactions between particles in those complicated moving interface problems, it could reasonably make physical significance if the Laplacian operator in (1) is replaced by its fractional version of Riesz-type potential to form the fractional Allen-Cahn equation. In this line, [16] proposed a kind of fractional Allen-Cahn model and discussed the solvability in some fractional Sobolev spaces, and [17] developed a fractional extension of the Allen-Cahn phase-field model with its fractional Laplacian defined by Riemann-Liouville fractional derivative that describes the mixture of two incompressible fluids. In [17], the authors also proposed a Petrov-Galerkin spectral method for spacial discretization combined with a stabilized ADI scheme for temporal discretization in the absence of rigorous numerical analysis for the solvability, stability, convergence, and conservation properties of the numerical scheme. Some other numerical methods such as finite difference [18], finite volume [19], and collocation method [20] were established. As far as we know, few research works have been done on the efficient finite element method and its rigorous numerical analysis for the fractional Laplace operator defined as Riesz-type potentials on the whole space R.In this article, we consider the following fractional Allen-Cahn equation [16]:(3)∂tu+K-Δru+ψ′u=0,x,t∈Ω×0,T,ux,0=u0x,x∈R,ux,t=0,x,t∈Ω¯×0,T,where Ω=(0,1), r∈(0,1),K>0 is the diffusion coefficient and the fractional Laplace operator (-Δ)r is defined as Riesz-type potentials on the whole space R in Section 2. The unknown u(x) can be viewed as an indicator of the concentration or volume fraction of one fluid at the location x in the immiscible mixture with the second fluid.The main objectives of this article are to (1) embed the fractional Allen-Cahn equation (1) into a Galerkin variational framework and thus develop its corresponding finite element procedure; (2) prove rigorously the solvability, the optimal convergence rates, and the physical properties, for examples, the mass decay, the energy dissipation, and the new energy equality, for the finite element solution; (3) combine the Toeplitz structure of the coefficient matrix and the merits of the classic CG algorithm [21] to design a fast CG (FCG) for the linearized finite element scheme, which reduces the computation cost and the storage to O(MlogM) and O(M), respectively; and (4) conduct numerical experiments to verify the efficiency of the FCG, which show that the FCG possesses the ideal convergence rates as Newtons algorithm does in space and time, preserves the mass, energy dissipation, and energy equality law, and recognizes accurately the phase separation by a very clear coarse graining process. The numerical experiments also test the tunable sharpness, that is, the influences of different fractional indices r and different coefficients η.The rest of this article is outlined as follows. Section2 is preliminaries. Section 3 is for the solvability and stability of the solution of discrete system. We demonstrate that the discrete solution preserves mass and energy dissipation and satisfies new energy equality in Section 4. Sections 5 and 6 are devoted to convergence analysis and efficient FCG algorithm, respectively. In the last section, numerical experiments are conducted to test the efficiency of the proposed efficient finite element algorithm.


## 2. Preliminaries

We first briefly revisit the definitions and some properties of fractional Laplace operator.Definition 1 (see [22, 23]).
Forr∈(0,1) and u∈S(x), where S(x) is the Schwartz class of rapidly decaying functions at infinity, the fractional Laplace operator (-Δ)r is reformulated by(4)-Δrux≔Crlimϵ→0∫R\Bx,ϵux-uyx-y1+2rdy,where Cr is a constant given by(5)Cr=4rrΓ1/2+rπ1/2Γ1-r.Definition 2 (see [22, 24, 25]).
For0<r<1, the fractional Sobolev spaces H0r(R) are defined by (6)L02R≔v∈L2R:v=0inR\Ω,H0rR≔v∈L02R:vx-vy2x-y1+2r∈L1R×R and equipped with the norm(7)vHrR2≔vL2R2+Cr2∫R∫Rvx-vy2x-y1+2rdxdy,v∈H0rR. and with equivalent seminorm(8)vHrR2≔Cr2∫R∫Rvx-vy2x-y1+2rdxdy,v∈H0rR. The energy space and the energy are defined, respectively, by(9)HrR≔v∈H0R:Erv<+∞,(10)Erv≔K2vH0r2+K4η2∫Rv4dx-K2η2∫Rv2dx,∀v∈H0rR.Definition 3 (see [22]).
For0<r<1, the operator Tr:H0r(R)→H0′r(R) is defined by(11)Tru,vH0r≔Cr2∫R∫Rux-uyvx-vyx-y1+2rdxdy,∀u,v∈H0r. Here H0′r(R) stands for its dual of H0r(R).It is easily seen thatTr is a symmetric positive definite operator.Definition 4 (see [22]).
LetT>0. If u∈C([0,T];H0r)∩W1,2([0,T];H0′r) satisfies(12)∂tu+K-Δru+1η2u3-u=0,a.e.inH0′r,then the u is called a weak solution to problem (3).


## 3. Finite Element Procedure

In this section, we construct finite element scheme for (3) based on the weak formulation (12) and prove the solvability and stability of solution of discrete system.Multiplying (12) by any v∈H0r, integrating over R and combining the homogeneous boundary condition, we obtain the variational formulation of (3) as to find u∈H0r such that(13)a∂tu,v+KTru,v+Kη2u3,v-Kη2u,v=0,∀v∈H0r,bux,0=u0x.TakingM and N as integers, we divide Ω=[0,1] uniformly by intervals Ii=[xi-1,xi] for i=1,2,…,M with x0=0, xM=1, h=hi=xi-xi-1=1/M, and partition the time interval [0,T] by the nodes tn=nτ for n=0,1,…,N with the time step τ=T/N.Upon the space partition, we define the finite element space as(14)Vh≔vh∈CΩ:vhIi∈PkIi,i=1,2,…,M;vh0=vh1=0;k≥1, here Pk(Ii) being the set of polynomials of degree not bigger than k over the interval Ii.Applying the backward Euler scheme to discrete the time derivative∂tu, we define the fully discrete finite element procedure of (13) as to find uhn∈Vh such that, for vh∈Vh,(15)auhn-uhn-1τ,vh+KTruhn,vh+Kη2uhn3,vh-Kη2uhn,vh=0,buh0=Rhu0,where Rhu0 is the elliptic projection of the initial value u0(x) to the finite element space.The dimension ofVh is (kM-1)×(kM-1); see the review [24]. Assume φi(x),i=1,2,…,kM-1 are the basis functions; then, we can express the numerical solution uhn by(16)uhnx=∑i=1kM-1uinφix,and, thus, the fully discrete finite element scheme (15) is transformed equivalently to the following algebraic equation:(17)1-τKη2A+τKBUn+τKη2HUnUn=AUn-1,where Un=(u1n,u2n,…,ukM-1n)T is the unknown vector, and the matrices A,B,H(Un) are(18)A=φi,φjkM-1×kM-1,B=Trφi,φjkM-1×kM-1,HUn=uhn2φi,φjkM-1×kM-1. It is easily verified that A, B, and H(Un) are symmetric and positive matrices.By using the contraction mapping principle, we prove the existence and uniqueness of the fully discrete finite element scheme (15).Theorem 5.
There exists a unique solutionUn,n=1,2,…,N to (17) for sufficiently small τ>0.Proof.
Selecting the time stepτ sufficiently small such that 0<τCK/η2<1 and noticing that τK>0 and the matrices A and B are positive definite, we know that the matrix (19)E≔1-τKη2A+τKB is positive definite and thus invertible for sufficiently small τ>0. Therefore, we solve Un from (17) (20)Un=-E-1τKη2HUnUn+E-1AUn-1.
Define the mappingT: X∈RkM-1→Y∈RkM-1 by(21)Y=-E-1τKη2HXX+E-1AUn-1. The mapping T is well defined for given Un-1 due to the positiveness of the matrix E.
We shall use a corollary of the well-known contraction mapping principle [26] to prove the mapping T having a unique fixed point in a bounded domain U={X∈RkM-1;X≤L} of RkM-1.
For this purpose, we letXi∈U and Yi=T(Xi) for i=1,2. Then,(22)Y1-Y2=TX1-TX2=-τKη2E-1HX1X1-HX2X2. Noticing that the matrix H(X) is Lipschitz continuously with respect to X, we obtain(23)Y1-Y2≤τKη2E-1HX1X1-HX2X2≤τKη2E-1HX1X1-X2+HX1-HX2X2≤τKη2E-1HX1X1-X2+CX1-X2X2=τKη2E-1HX1+CX2X1-X2≤τKη2E-1supX1,X2∈UHX1+CX2X1-X2. Select τ to be small enough such that(24)0<τ<Kη2E-1supX1,X2∈UHX1+CX2-1, and we have(25)Y1-Y2<X1-X2, which shows that the mapping T is a contractive mapping. This, together with an application of the well-known contraction mapping principle, completes the proof.Theorem 6.
Assume thatτ is small enough such that the conditions of Theorem 5 hold. Then, the fully discrete finite element scheme (15) is stable in the following sense, for J=1,2,…,N,(26)uhJ2+2τK∑i=1JuhiHr2+2τKη2∑i=1JuhiL44≤e2KT/η2uh02.Proof.
Takevh=uhn in (15) to obtain(27)auhn-uhn-1τ,uhn+KTruhn,uhn+Kη2uhn3,uhn-Kη2uhn,uhn=0,bu0h=Rhu0. Noting the fact that(28)Truhn,uhn=uhiHr2,uhn3,uhn=uhnL44 and(29)uhn-uhn-1,uhn≥12uhn2-12uhn-12, we have(30)uhn2-uhn-12+2τKuhiHr2+2τKη2uhnL44-uhn2≤0. Adding all the terms from n=1 to n=J, we get(31)uhJ2-uh02+2τK∑i=1JuhiHr2+2τKη2∑i=1JuhiL44-∑i=1Juhi2≤0. Namely,(32)uhJ2+2τK∑i=1JuhiHr2+2τKη2∑i=1JuhiL44≤2τKη2∑i=1Juhi2+uh02. Applying the discrete Gronwall inequality, we have(33)uhJ2+2τK∑i=1JuhiHr2+2τKη2∑i=1JuhiL44≤eJ2τK/η2uh02≤e2TK/η2uh02. This completes the proof.


## 4. Properties Preserved by the Finite Element Solution

In this section, we demonstrate that the finite element solution preserves the energy dissipation law and satisfies a redefined energy equality.Theorem 7.
The finite element solutionuhn of (15) preserves the energy dissipation law in the following sense, for J=1,2,…,N,(34)EruhJ≤EruhJ-1≤···≤Eru0.Proof.
Takingvh=uhn-uhn-1/τ in (15), we have(35)uhn-uhn-1τ,uhn-uhn-1τ+KTruhn,uhn-uhn-1τ+Kη2uhn3,uhn-uhn-1τ-Kη2uhn,uhn-uhn-1τ=0.Due to the fact that(36)uhnL22=uhn,uhn=UnTAUn,uhnHr2=Trun,un=UnTBUn, we have(37)uhn,uhn-uhn-1=UnTAUn-Un-1≥12UnTAUn-Un-1TAUn-1 and(38)Truhn,uhn-uhn-1=UnTBUn-Un-1≥12UnTBUn-Un-1TBUn-1. By applying the inequality ab≤1/2(a2+b2), we get(39)uhn3,uhn-uhn-1=uhn2,uhnuhn-uhn-1≥12uhn2,uhn2-uhn-12≥14uhn4-uhn-14,1. Then, combining the above inequalities with (35), we derive(40)τuhn-uhn-1τ2+K2uhnHr2-uhn-1Hr2+K4η2uhnL44-uhn-1L44-K2η2uhn2-uhn-12≤0.Therefore, we have(41)K2uhnHr2+K4η2uhnL44-K2η2uhn2≤K2uhn-1Hr2+K4η2uhn-1L44-K2η2uhn-12. It follows that Er(uhn)≤Er(uhn-1) for n=1,…,J. This completes the proof.We prove the following new energy equality preserved by the finite element solutionuh.Theorem 8.
The finite element solutionuhn of (15) satisfies the following energy equality for J=1,2,…,N,(42)EruhJ+1τ∑i=1Juhi-uhi-12+K2η2∑i=1Juhi-uhi-1Hr2+K4η2∑i=1Juhi2-uhi-122+K2η2∑i=1Juhi2-uhiuhi-12-K2η2∑i=1Juhi-uhi-12=Eruh0,where the left-hand side of the energy equality (42) is redefined as the new energy.Proof.
Choosinguhn-uhn-1/τ and multiplying τ2 to both sides in (15), we get(43)uhn-uhn-1,uhn-uhn-1+τKTruhn,uhn-uhn-1+τKη2uhn3,uhn-uhn-1-τKη2uhn,uhn-uhn-1=0.By a simple calculation, we have(44)Truhn,uhn-uhn-1=12uhnHr2-12uhn-1Hr2+12uhn-uhn-1Hr2,uhn,uhn-uhn-1=12uhn2-12uhn-12+12uhn-uhn-12 and(45)uhn3,uhn-uhn-1=14uhnL44-14uhn-1L44+14uhn2-uhn-122+12uhn2-uhnuhn-12. Substituting the above equalities into (43) and combining (10), we have(46)τEruhn-τEruhn-1+uhn-uhn-12+τK2η2uhn-uhn-1Hr2+τK4η2uhn2-uhn-122+τK2η2uhn2-uhnuhn-12-τK2η2uhn-uhn-12=0. Add all terms from n=1 to n=J to obtain (47)τEruhJ+∑i=1Juhi-uhi-12+τK2η2∑i=1Juhi-uhi-1Hr2+τK4η2∑i=1Juhi2-uhi-122+τK2η2∑i=1Juhi2-uhiuhi-12-τK2η2∑i=1Juhi-uhi-12=τEruh0. This completes the proof.


## 5. Convergence Analysis

In this section, we shall conduct convergence analysis for the fully discrete finite element scheme (15).Define the elliptic projectionRhu of the exact solution u as(48)Tru-Rhu,vh=0,∀vh∈Vh. The following estimate is valid for u∈Hs(Ω) and s>r [27],(49)u-Rhur≤Chmink+1-r,s-rus.Lett=tn and v=vh in (13), and subtract it from (15) to obtain the error equation(50)autn-uhn-uhn-1τ,vh+KTrun-uhn,vh+Kη2un3-uhn3,vh-Kη2un-uhn,vh=0,∀vh∈Vh,buh0=Rhu0.Denote (51)utn-un-un-1τ=Sn,un-Rhun=ρn,Rhun-uhn=θn, then (52)un-uhn=un-Rhun+Rhun-uhn=ρn+θn.Theorem 9.
Assume thatu∈H2(0,T;L2(R))∩H1(0,T;Hr(R)). Then, there exists integers m(1≤m≤N) and constants C0 and C1 independent of the step parameters h and τ such that(53)um-uhm2+2Kτ∑i=1mui-uhir2≤C0τ2+C1h2mink+1,s-r Here(54)C0=C0uttL20,T,L22,C1=C1utL20,T,Hs,uL20,T,Hs.Proof.
Takevh=θn in (50) to obtain(55)θn-θn-1τ,θn+KTrθn,θn+Kη2Rhun3-uhn3,θn=-Sn,θn+ρn-1-ρnτ,θn+Kη2Rhun3-un3,θn+Kη2ρn,θn+Kη2θn,θn.We first estimate the left-hand side of (55).
Applying theHo¨lder inequality and the Young inequality, we have(56)θn-θn-1τ,θn=1τθn,θn-θn-1,θn≥1τθn2-θn-122-θn22=12τθn2-θn-12. By using the monotonicity of function ((x1)3-(x2)3)(x1-x2)≥0, we have(57)Rhun3-uhn3,Rhun-uhn=Rhun3-uhn3,θn≥0. And then, we estimate the right-hand side of (55). Combining the Ho¨lder inequality and the Young inequality, we have(58)-Sn,θn≤-Snθn=1τ∫tn-1tnt-tn-1utttdtθn≤1τ∫01τ3∫tn-1tnutttdtdx1/2θn=τ1/2∫tn-1tnuttt2dt1/2θn≤τ∫tn-1tnuttt2dt+θn2. Similarly, we get(59)ρn-1-ρnτ,θn≤ρn-1-ρnτθn=1τ∫tn-1tnρttdtθn≤1τ1/2∫tn-1tnρtt2dt1/2θn≤12τ∫tn-1tnρtt2dt+12θn2. The third term on the right-hand side of (55) is estimated by(60)Rhun3-un3,θn≤C1ρnθn≤C12ρn2+C12θn2 and(61)ρn,θn≤ρnθn≤12ρn2+12θn2. Then, taking these inequalities into (55), we obtain(62)θn2-θn-12+2Kτθnr2≤τ2∫tn-1tnuttt2dt+∫tn-1tnρtt2dt+C1+1Kη2τρn2+2+C1+3Kη2τθn2. Noticing θ0=uh0-Rhu0=0 and adding all the terms from n=1 to n=N, we have(63)θm2+2Kτ∑i=1mθir2≤τ2∫0Tuttt2dt+∫0Tρtt2dt+C1K+Kη2τ∑i=1Nρi2+2η2+C1+3Kη2τ∑i=1Nθi2. The discrete Gronwall inequality is applied to (63), when τ is small enough to satisfy 2η2+C1+3K/η2τ<1,(64)θm2+2Kτ∑i=1mθir2≤e2η2+C1+3KT/η2τ2∫0Tuttt2dt+∫0Tρtt2dt+C1K+Kη2τ∑i=1Nρi2. If we take C=max{1/3e2η2+C1+3KT/η2,e2η2+C1+3KT/η2,C1K+K/η2e2η2+C1+3KT/η2}, then(65)θm2+2Kτ∑i=1mθir2≤Cτ2uttL20,T,L22+Ch2mink+1,sutL20,T,Hs2+Ch2mink+1,suL20,T,Hs2.
Then, using the elliptic projection estimate (49) and the triangle inequality, we obtain(66)um-uhm2+2Kτ∑i=1mui-uhir2≤ρm2+θm2+2Kτ∑i=1mθir2+2Kτ∑i=1mρir2≤C0τ2+C1h2mink+1,s-r.This concludes the proof.Remark 10.
Since our proof for the stability, energy dissipation, new energy-equality and error deduction of the discrete scheme (15) is independent of the dimensional argument, the analysis and conclusions in this article can be extended to multidimensional models, maybe with minor modification.


## 6. Fast-Conjugate-Gradient Algorithm (FCG)

The fully discrete finite element scheme (15) is nonlinear and can be solved by Newton’s iteration algorithm. Here we linearize the nonlinear terms by replacing (uhn)2 by (uhn-1)2 and solve this linearized version by the conjugated gradient algorithm (CG). We find that if the CG algorithm is directly used, the computation cost and storage will reach up to O(M2) since the coefficient matrix B is nonsparse due to the nonlocality of the fractional Laplace operator. We also find that if the linear finite element space is employed, the matrix B is a Toeplitz matrix (see Section 7.1), which makes it possible to reduce the computation cost and storage to O(Mlog⁡M) by a delicate combination of the CG algorithm, the fast Fourier transform (FFT) and its Toeplitz structure of the matrix B.In this section, we depict the general ideas for the combination of CG, FFT, and the Toeplitz matrix and design a fast CG algorithm (FCG). For a good review of a fast algorithm generating from FFT and the Toeplitz matrix, we refer to [28].Noticing what causes the computation cost up toO(M2) is a direct use of the matrix-vector multiplication in the classical CG algorithm; as a remedy, we can obtain the matrix-vector multiplication BkM-1x for a vector x∈RkM-1 and (kM-1)×(kM-1) Toeplitz matrix BkM-1 through the following two steps [29, 30]:(1) assembling the(kM-1)×(kM-1) Toeplitz matrix BkM-1 into a 2(kM-1)×2(kM-1) cyclic matrix C2(kM-1);(2) extracting the Toeplitz matrix-vector multiplicationBkM-1x for a vector x∈RkM-1 via(67)C2kM-1X=BkM-1DkM-1DkM-1BkM-1x0=BkM-1xDkM-1x.In (67), DkM-1 is determined simply by the matrix BkM-1 and the computational cost and the storage of the matrix-vector multiplication BkM-1x are O(MlogM) and O(M) [30]. This will reduce the computation cost and the storage from O(M2) to O(MlogM) and O(M) compared to the direct use of matrix-vector multiplication.Upon this analysis, we write down the improved algorithm, called the fast CG (FCG), sentence by sentence in Algorithm1, and present the corresponding conclusion by Theorem 11.Algorithm 1:The FCG Algorithm.
function FCG(F,U)           w=FT(p); rb=F-FT(U);            α=ρ/(w′∗p); nr=norm(rb);              U=U+α∗p; k=0;                rb=rb-α∗w; ρ=rb′∗rb;               ρwave=ρ; while1                ρ=rb′∗rb; k=k+1;              if  sqrt(rho)>e-7  then if  k==1  then                continue; p=rb;                else else                    break; β=ρ/ρwave;              end if p=rb+β∗p;         end function end ifTheorem 11.
The computational cost and storage of the fast conjugated gradient algorithm (FCG) are reduced fromO(M2) and O(M2) to O(MlogM) and O(M) compared with the traditional CG algorithm.


## 7. Numerical Experiments

In this section, we linearize the nonlinear term in the fully discrete finite element scheme (15) by taking (uhn-1)2uhn to replace (uhn)3; we carry out two numerical experiments to verify the result convergence analysis and the physical property of the solution of discrete system.


### 7.1. The Linear Finite Element System

For simplicity, we take the linear finite element system to compute algebraic equation (17).Fork=1, the inner nodal basis function is φi(x) at xi,i=1,2,…,M-1, with the following structure:(68)φix=x-xi-1h,x∈xi-1,xi,xi+1-xh,x∈xi,xi+1,0,else. For i,j=1,2,…,M-1, the matrix A is(69)A=φi,φjM-1×M-1=2h3h60⋯0h62h3h6⋯0⋮⋱⋱⋱⋮0⋱⋱2h3h600⋯h62h3. We replace the nonlinear term H(Un) with H(Un-1),(70)HUn-1=uhn-12φi,φjM-1×M-1,i,j=1,2,…,M-1. After calculation, the ij-entries of the matrix H(Un-1) are(71)uhn-12φi,φj=ui-1n-1230+2uin-125+ui-1n-1uin-1+uin-1ui+1n-110+ui+1n-1230,j-i=0,uin-1220+ujn-1220+uin-1ujn-115,j-i=1,0,j-i=m,m=2,…,M-2.For i,j=1,2,…,M-1, the entries of the matrix B are respectively(72)Trφi,φj=Cr2∫R∫Rφix-φiyφjx-φjyx-y1+2rdxdy=Cr2∫Ω∫R\Ωφix-φiyφjx-φjyx-y1+2rdxdy+∫R\Ω∫Ωφix-φiyφjx-φjyx-y1+2rdxdy+∫Ω∫Ωφix-φiyφjx-φjyx-y1+2rdxdy.After calculating the above integrals (72), we obtain that the ij-entries of matrix B are solved for r≠1/2(73)Trφi,φj=h1-2rD-4+23-2r,j-i=0,h1-2rD72-24-2r+33-2r2,j-i=1,h1-2rD3m3-2r-2m+13-2r-2m-13-2r+m-23-2r2+m+23-2r2,j-i=m,m=2,…,M-2,where (74)D=r1-2r2-2r3-2rCr,Cr=4rrΓ1/2+rπ1/2Γ1-r. For r=1/2,(75)Trφi,φj=4Crln2,j-i=0,Cr72ln3-ln2,j-i=1,Cr14ln2-9ln3,j-i=2,Cr3m2lnm+m-222lnm-2+m+222lnm+2-2m-12lnm-1-2m+12lnm+1,j-i=m,m=3,…,M-2.


### 7.2. Tests on the Efficiency of the Finite Element Procedure and the FCG

Example 12.
AssumingΩ=[0,1],T=1,K=1, the analytic solution is(76)ux,t=x1-xet∈H01+γΩ, where γ∈(0,1/2), is selected as close to 1/2 as possible. The source term f should be expressed as follows:(77)f=x1-xet+KCrx1-2ret2r1-2r+1-x1-2ret2r1-2r-x2-2retr1-2r2-2r-1-x2-2retr1-2r2-2r+Kη2x31-x3et-x1-xet.
We apply Example12 to verify the error convergence rate of the linear finite element system and to test the efficiency of the FCG for r=1/3,2/3.
The numerical solutionUh(x) is calculated by the FCG for T=1. The numerical solution and the exact solution are identical at t=1 in Figure 1, which show that the FCG is accurate.
Table1 show the spatial convergence rates with respect to u-uhH0r for r=1/3,2/3 respectively. By comparing the convergence rates of FCG to the convergence rate of Newton, Gauss, and CG, respectively, we conclude that the FCG algorithm is efficient. The convergence rates of u-uhH0r are at least 1.6 and 1.2 for r=1/3,2/3, respectively, which are almost equal to the theoretical error rate min{1+γ,2}-r.
Table2 shows the spatial convergence rates with respect to u-uh for r=1/3,2/3, respectively. The convergence rates of u-uh are at least 1.9 and 1.5 for r=1/3,2/3, respectively. We observe that the spatial convergence rates decrease as the spatial fractional order r increases.Table 1
Spatial errors and convergence rates foru-uhHr,0.
r h Newton rate Gauss rate CG rate FCG rate 1 3 2 - 3 1.247E-2  1.336E-2  1.336E-2  1.336E-2 2 - 4 4.026E-3 1.63 4.204E-3 1.66 4.204E-3 1.66 4.204E-3 1.66 2 - 5 1.276E-3 1.65 1.312E-3 1.67 1.312E-3 1.67 1.312E-3 1.67 2 - 6 4.047E-4 1.65 4.121E-4 1.67 4.121E-4 1.67 4.121E-4 1.67 2 - 7 1.299E-4 1.63 1.313E-4 1.64 1.313E-4 1.64 1.313E-4 1.64   2 3 2 - 3 4.260E-2  4.294E-2  4.294E-2  4.294E-2 2 - 4 1.808E-2 1.23 1.814E-2 1.24 1.814E-2 1.24 1.814E-2 1.24 2 - 5 7.544E-3 1.26 7.556E-3 1.26 7.556E-3 1.26 7.556E-3 1.26 2 - 6 3.164E-3 1.25 3.167E-3 1.25 3.167E-3 1.25 3.167E-3 1.25 2 - 7 1.361E-3 1.21 1.362E-3 1.21 1.362E-3 1.21 1.362E-3 1.21Table 2
Spatial errors and convergence rates foru-uh.
r h Newton rate Gauss rate CG rate FCG rate 1 3 2 - 3 6.117E-3  7.269E-3  7.270E-3  7.434E-3 2 - 4 1.755E-3 1.80 2.031E-3 1.84 2.030E-3 1.84 2.060E-3 1.85 2 - 5 4.800E-4 1.87 5.476E-4 1.89 5.475E-4 1.89 5.540E-4 1.89 2 - 6 1.286E-4 1.90 1.451E-4 1.92 1.451E-4 1.92 1.468E-4 1.92 2 - 7 3.434E-5 1.90 3.834E-5 1.92 3.833E-5 1.92 3.891E-5 1.92   2 3 2 - 3 6.517E-3  7.322E-3  7.321E-3  7.497E-3 2 - 4 2.220E-3 1.55 2.411E-3 1.60 2.414E-3 1.60 2.452E-3 1.61 2 - 5 7.480E-4 1.56 7.937E-4 1.60 7.937E-4 1.60 8.050E-4 1.60 2 - 6 2.567E-4 1.54 2.675E-4 1.56 2.675E-4 1.56 2.711E-4 1.57 2 - 7 9.084E-5 1.49 9.340E-5 1.51 9.340E-5 1.51 9.455E-5 1.51Figure 1
Comparison of the initial value, the exact solution, and the numerical solution forr=1/3,K=1,η=1.Remark 13.
RecallingL2 error estimate for second order elliptic equation, one can obtain the estimate of u-uh one order higher than that for u-uh1. However, in the case of the fractional order elliptic equation, we only obtain the convergence rate of u-uh dependent on the index r. From the numerical Table 2, the convergence order of u-uh conforms to the following relation.Figure2 shows that the convergence rate of u-uh dependent on the index r, as the index r increases, the rates of spatial convergence decreases and then increases.Figure 2
The relation of spatial convergence rates and fractional orderr.Table3 shows the temporal convergence rates as r=1/3,r=2/3, respectively. The temporal convergence rates are almost equal to 1 for r=1/3,r=2/3 at h=2-10, which are consistent with the theoretical expectation for the backward Euler scheme.Table 3
Temporal errors and convergence rates.
r τ Newton rate Gauss rate CG rate FCG rate 1 3 1 8 1.173E-2  2.245E-2  2.245E-2  2.246E-2 1 10 6.012E-3 0.96 1.788E-2 1.020 1.788E-2 1.020 1.788E-2 1.020 1 12 3.044E-3 0.98 1.486E-2 1.017 1.486E-2 1.017 1.486E-2 1.017 1 14 1.532E-3 0.99 1.270E-2 1.014 1.270E-2 1.014 1.270E-2 1.014 1 16 7.696E-4 0.99 1.110E-2 1.012 1.110E-2 1.012 1.110E-2 1.012   2 3 1 8 6.768E-3  1.375E-2  1.376E-2  1.376E-2 1 10 3.478E-3 0.96 1.103E-2 0.988 1.103E-2 0.989 1.103E-2 0.988 1 12 1.766E-3 0.97 9.210E-3 0.991 9.211E-3 0.991 9.211E-3 0.991 1 14 8.935E-4 0.98 7.904E-3 0.992 7.904E-3 0.992 7.904E-3 0.992 1 16 4,527E-4 0.98 6.922E-3 0.993 6.922E-3 0.993 6.922E-3 0.993Table4 shows the efficiency of the FCG algorithm. It is easily obtained that the CPU time consumed by the FCG is much less than that of the Newton, Gauss, and CG as h decreases. It can be seen that the CPU time of FCG is 25.35 minutes while the Newton’s 15.2 hours and the CG’s 6.79 hours for h=2-10, r=1/3. These show that the FCG is an efficiency algorithm.Table 4
CPU time of Newton, Gauss, CG, and FCG.
h Newton Gauss CG FCG r=1/3 2 - 5 0.38s 0.25s 0.22s 0.34s 2 - 6 5.11s 8.22s 2.14s 1.12s 2 - 7 1.1m 3.77m 16.63s 4.94s 2 - 8 19.3m 1.94h 2.27m 14.71s 2 - 9 5.19h 69.9h 23.6m 1.74m 2 - 10 15.2h - 6.79h 25.35m


### 7.3. Tests on Physical Properties of Finite Element Solution

Example 14.
Lettingh=1/100, τ=1/500, the initial value is(78)ux,0=0.9sin2πx.
This example is used to test the physical properties such as the energy laws and the new energy equality, the mass dissipation, the phase separation, and the tunable sharpness in Figures3–8.
We present the energy dissipation in Figure3 and the new energy equality in Figure 4 for r=0.3, η=0.85, and K=1. In Figure 3, the energy dissipates to zero as time increases. In Figure 4, the energy jumps in scale 10-7 at the beginning until t=1 and keeps the same at the end; this demonstrates that the finite element solution to (15) almost satisfies new equality equation.
We plot the mass changes forr=0.3, η=0.85, and K=1 in Figure 5. The mass value dissipates as time increases; this property is coinciding with the classical Allen-Cahn equation, and this also explains the applicability of fractional Laplacian for Allen-Cahn equation in some way.
We present the physical phenomenon of a coarse graining process as the time increases forr=0.6, η=0.1, and K=0.001 in Figure 6. The two phases are connected in a smooth steep line; this indicates that the phase separation is complete.
Figure7 shows that the interface gets sharper and sharper as r decreases. Figure 8 shows that the interface gets sharper and sharper as η decreases. The above two groups of experimental results show that the tunable sharpness of the interface is related to different parameters r and η.Figure 3
Energy dissipation with time increasing.Figure 4
The new energy equality forr=0.3, η=0.85, and K=1.Figure 5
The mass dissipation forr=0.3, η=0.85, and K=1.Figure 6
The coarse graining process with respect to differentt.Figure 7
The tunable sharpness with respect to differentr.Figure 8
The tunable sharpness with respect to differentη.


## 7.1. The Linear Finite Element System

For simplicity, we take the linear finite element system to compute algebraic equation (17).Fork=1, the inner nodal basis function is φi(x) at xi,i=1,2,…,M-1, with the following structure:(68)φix=x-xi-1h,x∈xi-1,xi,xi+1-xh,x∈xi,xi+1,0,else. For i,j=1,2,…,M-1, the matrix A is(69)A=φi,φjM-1×M-1=2h3h60⋯0h62h3h6⋯0⋮⋱⋱⋱⋮0⋱⋱2h3h600⋯h62h3. We replace the nonlinear term H(Un) with H(Un-1),(70)HUn-1=uhn-12φi,φjM-1×M-1,i,j=1,2,…,M-1. After calculation, the ij-entries of the matrix H(Un-1) are(71)uhn-12φi,φj=ui-1n-1230+2uin-125+ui-1n-1uin-1+uin-1ui+1n-110+ui+1n-1230,j-i=0,uin-1220+ujn-1220+uin-1ujn-115,j-i=1,0,j-i=m,m=2,…,M-2.For i,j=1,2,…,M-1, the entries of the matrix B are respectively(72)Trφi,φj=Cr2∫R∫Rφix-φiyφjx-φjyx-y1+2rdxdy=Cr2∫Ω∫R\Ωφix-φiyφjx-φjyx-y1+2rdxdy+∫R\Ω∫Ωφix-φiyφjx-φjyx-y1+2rdxdy+∫Ω∫Ωφix-φiyφjx-φjyx-y1+2rdxdy.After calculating the above integrals (72), we obtain that the ij-entries of matrix B are solved for r≠1/2(73)Trφi,φj=h1-2rD-4+23-2r,j-i=0,h1-2rD72-24-2r+33-2r2,j-i=1,h1-2rD3m3-2r-2m+13-2r-2m-13-2r+m-23-2r2+m+23-2r2,j-i=m,m=2,…,M-2,where (74)D=r1-2r2-2r3-2rCr,Cr=4rrΓ1/2+rπ1/2Γ1-r. For r=1/2,(75)Trφi,φj=4Crln2,j-i=0,Cr72ln3-ln2,j-i=1,Cr14ln2-9ln3,j-i=2,Cr3m2lnm+m-222lnm-2+m+222lnm+2-2m-12lnm-1-2m+12lnm+1,j-i=m,m=3,…,M-2.


## 7.2. Tests on the Efficiency of the Finite Element Procedure and the FCG

Example 12.
AssumingΩ=[0,1],T=1,K=1, the analytic solution is(76)ux,t=x1-xet∈H01+γΩ, where γ∈(0,1/2), is selected as close to 1/2 as possible. The source term f should be expressed as follows:(77)f=x1-xet+KCrx1-2ret2r1-2r+1-x1-2ret2r1-2r-x2-2retr1-2r2-2r-1-x2-2retr1-2r2-2r+Kη2x31-x3et-x1-xet.
We apply Example12 to verify the error convergence rate of the linear finite element system and to test the efficiency of the FCG for r=1/3,2/3.
The numerical solutionUh(x) is calculated by the FCG for T=1. The numerical solution and the exact solution are identical at t=1 in Figure 1, which show that the FCG is accurate.
Table1 show the spatial convergence rates with respect to u-uhH0r for r=1/3,2/3 respectively. By comparing the convergence rates of FCG to the convergence rate of Newton, Gauss, and CG, respectively, we conclude that the FCG algorithm is efficient. The convergence rates of u-uhH0r are at least 1.6 and 1.2 for r=1/3,2/3, respectively, which are almost equal to the theoretical error rate min{1+γ,2}-r.
Table2 shows the spatial convergence rates with respect to u-uh for r=1/3,2/3, respectively. The convergence rates of u-uh are at least 1.9 and 1.5 for r=1/3,2/3, respectively. We observe that the spatial convergence rates decrease as the spatial fractional order r increases.Table 1
Spatial errors and convergence rates foru-uhHr,0.
r h Newton rate Gauss rate CG rate FCG rate 1 3 2 - 3 1.247E-2  1.336E-2  1.336E-2  1.336E-2 2 - 4 4.026E-3 1.63 4.204E-3 1.66 4.204E-3 1.66 4.204E-3 1.66 2 - 5 1.276E-3 1.65 1.312E-3 1.67 1.312E-3 1.67 1.312E-3 1.67 2 - 6 4.047E-4 1.65 4.121E-4 1.67 4.121E-4 1.67 4.121E-4 1.67 2 - 7 1.299E-4 1.63 1.313E-4 1.64 1.313E-4 1.64 1.313E-4 1.64   2 3 2 - 3 4.260E-2  4.294E-2  4.294E-2  4.294E-2 2 - 4 1.808E-2 1.23 1.814E-2 1.24 1.814E-2 1.24 1.814E-2 1.24 2 - 5 7.544E-3 1.26 7.556E-3 1.26 7.556E-3 1.26 7.556E-3 1.26 2 - 6 3.164E-3 1.25 3.167E-3 1.25 3.167E-3 1.25 3.167E-3 1.25 2 - 7 1.361E-3 1.21 1.362E-3 1.21 1.362E-3 1.21 1.362E-3 1.21Table 2
Spatial errors and convergence rates foru-uh.
r h Newton rate Gauss rate CG rate FCG rate 1 3 2 - 3 6.117E-3  7.269E-3  7.270E-3  7.434E-3 2 - 4 1.755E-3 1.80 2.031E-3 1.84 2.030E-3 1.84 2.060E-3 1.85 2 - 5 4.800E-4 1.87 5.476E-4 1.89 5.475E-4 1.89 5.540E-4 1.89 2 - 6 1.286E-4 1.90 1.451E-4 1.92 1.451E-4 1.92 1.468E-4 1.92 2 - 7 3.434E-5 1.90 3.834E-5 1.92 3.833E-5 1.92 3.891E-5 1.92   2 3 2 - 3 6.517E-3  7.322E-3  7.321E-3  7.497E-3 2 - 4 2.220E-3 1.55 2.411E-3 1.60 2.414E-3 1.60 2.452E-3 1.61 2 - 5 7.480E-4 1.56 7.937E-4 1.60 7.937E-4 1.60 8.050E-4 1.60 2 - 6 2.567E-4 1.54 2.675E-4 1.56 2.675E-4 1.56 2.711E-4 1.57 2 - 7 9.084E-5 1.49 9.340E-5 1.51 9.340E-5 1.51 9.455E-5 1.51Figure 1
Comparison of the initial value, the exact solution, and the numerical solution forr=1/3,K=1,η=1.Remark 13.
RecallingL2 error estimate for second order elliptic equation, one can obtain the estimate of u-uh one order higher than that for u-uh1. However, in the case of the fractional order elliptic equation, we only obtain the convergence rate of u-uh dependent on the index r. From the numerical Table 2, the convergence order of u-uh conforms to the following relation.Figure2 shows that the convergence rate of u-uh dependent on the index r, as the index r increases, the rates of spatial convergence decreases and then increases.Figure 2
The relation of spatial convergence rates and fractional orderr.Table3 shows the temporal convergence rates as r=1/3,r=2/3, respectively. The temporal convergence rates are almost equal to 1 for r=1/3,r=2/3 at h=2-10, which are consistent with the theoretical expectation for the backward Euler scheme.Table 3
Temporal errors and convergence rates.
r τ Newton rate Gauss rate CG rate FCG rate 1 3 1 8 1.173E-2  2.245E-2  2.245E-2  2.246E-2 1 10 6.012E-3 0.96 1.788E-2 1.020 1.788E-2 1.020 1.788E-2 1.020 1 12 3.044E-3 0.98 1.486E-2 1.017 1.486E-2 1.017 1.486E-2 1.017 1 14 1.532E-3 0.99 1.270E-2 1.014 1.270E-2 1.014 1.270E-2 1.014 1 16 7.696E-4 0.99 1.110E-2 1.012 1.110E-2 1.012 1.110E-2 1.012   2 3 1 8 6.768E-3  1.375E-2  1.376E-2  1.376E-2 1 10 3.478E-3 0.96 1.103E-2 0.988 1.103E-2 0.989 1.103E-2 0.988 1 12 1.766E-3 0.97 9.210E-3 0.991 9.211E-3 0.991 9.211E-3 0.991 1 14 8.935E-4 0.98 7.904E-3 0.992 7.904E-3 0.992 7.904E-3 0.992 1 16 4,527E-4 0.98 6.922E-3 0.993 6.922E-3 0.993 6.922E-3 0.993Table4 shows the efficiency of the FCG algorithm. It is easily obtained that the CPU time consumed by the FCG is much less than that of the Newton, Gauss, and CG as h decreases. It can be seen that the CPU time of FCG is 25.35 minutes while the Newton’s 15.2 hours and the CG’s 6.79 hours for h=2-10, r=1/3. These show that the FCG is an efficiency algorithm.Table 4
CPU time of Newton, Gauss, CG, and FCG.
h Newton Gauss CG FCG r=1/3 2 - 5 0.38s 0.25s 0.22s 0.34s 2 - 6 5.11s 8.22s 2.14s 1.12s 2 - 7 1.1m 3.77m 16.63s 4.94s 2 - 8 19.3m 1.94h 2.27m 14.71s 2 - 9 5.19h 69.9h 23.6m 1.74m 2 - 10 15.2h - 6.79h 25.35m


## 7.3. Tests on Physical Properties of Finite Element Solution

Example 14.
Lettingh=1/100, τ=1/500, the initial value is(78)ux,0=0.9sin2πx.
This example is used to test the physical properties such as the energy laws and the new energy equality, the mass dissipation, the phase separation, and the tunable sharpness in Figures3–8.
We present the energy dissipation in Figure3 and the new energy equality in Figure 4 for r=0.3, η=0.85, and K=1. In Figure 3, the energy dissipates to zero as time increases. In Figure 4, the energy jumps in scale 10-7 at the beginning until t=1 and keeps the same at the end; this demonstrates that the finite element solution to (15) almost satisfies new equality equation.
We plot the mass changes forr=0.3, η=0.85, and K=1 in Figure 5. The mass value dissipates as time increases; this property is coinciding with the classical Allen-Cahn equation, and this also explains the applicability of fractional Laplacian for Allen-Cahn equation in some way.
We present the physical phenomenon of a coarse graining process as the time increases forr=0.6, η=0.1, and K=0.001 in Figure 6. The two phases are connected in a smooth steep line; this indicates that the phase separation is complete.
Figure7 shows that the interface gets sharper and sharper as r decreases. Figure 8 shows that the interface gets sharper and sharper as η decreases. The above two groups of experimental results show that the tunable sharpness of the interface is related to different parameters r and η.Figure 3
Energy dissipation with time increasing.Figure 4
The new energy equality forr=0.3, η=0.85, and K=1.Figure 5
The mass dissipation forr=0.3, η=0.85, and K=1.Figure 6
The coarse graining process with respect to differentt.Figure 7
The tunable sharpness with respect to differentr.Figure 8
The tunable sharpness with respect to differentη.



---

*Source: 7969371-2019-07-28.xml*