# A Simple Benchmark Problem for the Numerical Methods of the Cahn–Hilliard Equation

Authors: Yibao Li; Chaeyoung Lee; Jian Wang; Sungha Yoon; Jintae Park; Junseok Kim
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

We present a very simple benchmark problem for the numerical methods of the Cahn–Hilliard (CH) equation. For the benchmark problem, we consider a cosine function as the initial condition. The periodic sinusoidal profile satisfies both the homogeneous and periodic boundary conditions. The strength of the proposed problem is that it is simpler than the previous works. For the benchmark numerical solution of the CH equation, we use a fourth-order Runge–Kutta method (RK4) for the temporal integration and a centered finite difference scheme for the spatial differential operator. Using the proposed benchmark problem solution, we perform the convergence tests for an unconditionally gradient stable scheme via linear convex splitting proposed by Eyre and the Crank–Nicolson scheme. We obtain the expected convergence rates in time for the numerical schemes for the one-, two-, and three-dimensional CH equations.


---

## Body

## 1. Introduction

We present a very simple benchmark problem for the numerical methods of the following Cahn–Hilliard (CH) equation [1–3]:(1)∂ϕx,t∂t=ΔF′ϕx,t−ε2Δϕx,t,x,t∈Ω×0,T,where Ω⊂ℝdd=1,2,3 is a bounded domain, ϕx,t is a compositional field, Fϕ=0.25ϕ2−12, and ε is a positive constant. The CH equation was proposed for a model of the spinodal decomposition in a binary mixture and has been used to model many scientific phenomena such as topology optimization [4], phase separation [5–7], image processing [8], two-phase fluid flows [9, 10], crystal model [11], tumor growth [12, 13], and microstructure formations (see [14] for the basic principles and practical applications and [3] for the physical, mathematical, and numerical derivations of the CH equation).There have been many research papers regarding the accurate numerical methods for the CH equation [15–18] (see the references therein for more details). In particular, a comparison study among various numerical methods for the CH equation was conducted on the perspective of stability and efficiency [19]. A test set developed to verify models, algorithms, or the accuracy of numerical methods is called “benchmark problem” [20]. However, there are only few benchmark problems for validating the accuracy of the proposed numerical methods. Recently, the authors in [21] proposed a benchmark problem for the two- and three-dimensional CH equations. The benchmark problems were the shrinking annulus and spherical shell in the two- and three-dimensional cases, respectively. They used the explicit Euler scheme with a very fine time step size, Δt<10−5. The authors in [22] presented four benchmark problems for the Allen–Cahn (AC) and CH equations. The benchmark is the time T at which the value at a point in the domain changes from negative to positive. The authors in [20] presented two benchmark problems for phase-field models of solute diffusion and phase separation. Recently, a verification method for the convergence rates of the numerical solutions for well-known parabolic partial differential equations was proposed in [23].Typically, benchmark testing is essential for comparing the performance of a method. However, many researchers conduct self-tests using their own methods in most cases. There is a problem such that an error may occur relatively compared to the actual solution because an approximated solution converges to a numerical solution via an employed numerical method although convergence of a solution through numerical analysis is guaranteed. Moreover, most of these testing schemes are often implicit methods, which have relatively low accuracy compared to the results of long time simulation of explicit methods with a small time step. Therefore, the main purpose of this paper is to present a very simple benchmark problem for the numerical methods of the CH equation, which does not employ the self-test but a classical explicit method for the temporal discretization. The strength of the proposed problem is that it is simpler than the previous works [20–22]. Furthermore, the benchmark problem can be applied to both Neumann and periodic boundary conditions.The contents of this paper are summarized as follows. In Section2, a numerical solution is presented. In Section 3, numerical results are shown. Conclusions are made in Section 4.


## 2. Numerical Solution

We consider two discretization schemes to the CH equation for benchmark problem; one is the linearly stabilized splitting scheme (LSS) proposed by Eyre [24] which is first-order accurate in time and the other is the Crank–Nicolson (CN) scheme [25] which is second-order in time. In order to perform a benchmark test for these two methods, we first need to obtain a benchmark numerical solution. To obtain the benchmark solution, we consider the one-dimensional CH equation:(2)ϕtx,t=F′ϕx,t−ε2ϕxxx,txx,x,t∈0,1×0,T,with the homogeneous Neumann boundary condition ϕx0,t=ϕx1,t=ϕxxx0,t=ϕxxx1,t=0. Let h=1/Nx be the uniform spatial step size, xi=i−0.5h for 1≤i≤Nx, and tn=nΔtref, where Δtref=T/Nt is the reference time step and T is the final time. We define ϕin as an approximation of ϕxi,tn. Now, we discretize equation (2) in time using a fourth-order Runge–Kutta method (RK4) [26] as follows:(3)αin=ΔtrefΔhF′ϕin−ε2Δhϕin,βin=ΔtrefΔhF′ϕin+0.5αin−ε2Δhϕin+0.5αin,γin=ΔtrefΔhF′ϕin+0.5βin−ε2Δhϕin+0.5βin,δin=ΔtrefΔhF′ϕin+γin−ε2Δhϕin+γin,ϕin+1=ϕin+16αin+2βin+2γin+δin,fori=1,2,…,Nx,where Δhψi=ψi−1−2ψi+ψi+1/h2 is the standard discrete Laplacian with ψ0=ψ1 and ψNx+1=ψNx.


## 3. Numerical Results

### 3.1. Convergence and Stability Tests

First, we conduct a traditional convergence test. As a concrete example, an initial condition is given by(4)ϕx,0=0.1cos2πxonx∈Ω=0,1.Note that equation (4) can be used in testing numerical schemes with both the Neumann and periodic boundary conditions. We employ ε=16h/22tanh−10.9≈0.03, which is the approximate transition layer width [27]. Table 1 shows the rate of convergence. We measure the accuracy of time at the final time T using eΔt,0.5Δt2 which is defined as follows:(5)eΔt,Δt/22=1Nx∑i=1NxeiΔt,Δt/22,where eiΔt,Δt/2 is an error between the computed solutions using large and small time steps, respectively. More precisely, it is(6)eiΔt,Δt/2=ϕiΔt−ϕiΔt/2.Table 1
Rate of convergence of RK4.
Δt,Δt/2h4,h4/2Rateh4/2,h4/4Rateh4/4,h4/8eΔt,Δt/224.378e–114.0302.681e–124.0971.567e–13In this test, we adopt the space step size ash=1/10 and the final time as T=10h4. We can confirm that RK4 is of fourth-order in time indeed according to Table 1.For the next test, we investigate the maximum time step sizesΔtc with respect to the spatial step size h. Here, we fix Nt=105 and other parameters and the initial condition are the same as those described above. Table 2 lists the largest time step sizes satisfying that the numerical solution does not blow up.Table 2
Practical maximum time step sizeΔtc with respect to h.
h1/161/321/641/128Δtc4.686e−51.152e−52.810e−67.085e−7Figure1 depicts the practically stable region and maximum time step size Δtc as described in Table 2. The explicit scheme is theoretically Δt≈Oh4.Figure 1
Practical maximum time step sizeΔtc with respect to h.Figure2 represents the CPU times (seconds) with respect to the number of spatial step Nx, performed on Intel Core i5-6400 CPU at 2.70 GHz with 4 GB RAM. When we vary Nx as 16, 32, 64, and 128 with fixed parameters, Δt=1/1284 and Nt=107, the CPU times are 88.812, 175.468, 354.14, and 708.718. That is, when Nx is doubled, the CPU time is approximately doubled.Figure 2
CPU times versus spatial step.


### 3.2. Benchmark Problem

Let us consider the benchmark problem for the CH equation. We adopt equation (4) as an initial condition. For the reference numerical solution of the benchmark problem, we use h=1/128, Δtref=h4, and Nt=107. Therefore, the final time is T=NtΔtref=107h4. For the model parameter, we use ε=24h/22tanh−10.9≈0.045. Figure 3 shows the initial condition and the reference numerical solution at t=107h4.Figure 3
Initial condition and the reference numerical solution att=107h4.


#### 3.2.1. One-Dimensional Space

Let us define an errorex,T as the difference between ϕΔt and ϕΔtref at the final time T, i.e., ex,T=ϕΔtx,T−ϕΔtrefx,T, where ϕΔt and ϕΔtref are the numerical solution with a time step Δt and the reference numerical solution with a time step Δtref, respectively. The one-dimensional LSS and CN schemes to be benchmarked are as follows:(7)ϕin+1−ϕinΔt=Δhϕin3−3ϕin+2ϕin+1−ε2Δhϕin+1,(8)ϕin+1−ϕinΔt=12Δhϕin3−ϕin−ε2Δhϕin+ϕin+13−ϕin+1−ε2Δhϕin+1,where Δt is a time step. We solve equations (7) and (8) using a multigrid method [25, 28]. Note that we use the Gauss–Seidel iterative method in multigrid process combined with Newton’s approximation to compute the nonlinear term in equation (8). Figures 4(a) and 4(b) show the errors for the LSS and CN scheme, respectively.Figure 4
Errors between the reference solution and the solutions evaluated using (a) the LSS and (b) the CN scheme att=107h4.
(a)(b)Now, we define the discretel2-norm and the discrete maximum norm as(9)e2=1Nx∑i=1Nxei2,e∞=max1≤i≤Nxei,respectively. Tables 3 and 4 list the numerical errors and convergence rates of the LSS and those of the CN scheme at t=107h4, respectively. We confirm that both the schemes achieve the expected convergence rates (the first-order for the LSS and the second-order for the CN scheme).Table 3
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt1250h4Rate625h4Rate312.5h4Rate156.25h4e25.687e − 40.9952.854e − 40.9971.430e − 40.9997.150e − 5e∞8.970e − 40.9954.500e − 40.9982.254e − 40.9991.128e − 4Table 4
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 92.0011.600e − 9e∞1.608e − 72.0004.020e − 82.0001.010e − 82.0002.500e − 9


#### 3.2.2. Two-Dimensional Space

Next, we consider the two-dimensional version of equations (7) and (8) to the CH equation on Ω=0,1×0,1. Straightforward extensions of the LSS and CN schemes are as follows:(10)ϕijn+1−ϕijnΔt=Δhϕijn3−3ϕijn+2ϕijn+1−ε2Δhϕijn+1,(11)ϕijn+1−ϕijnΔt=12Δhϕijn3−ϕijn−ε2Δhϕijn+ϕijn+13−ϕijn+1−ε2Δhϕijn+1,where Δhψij=ψi−1,j+ψi+1,j+ψi,j−1+ψi,j+1−4ψij/h2 and h=1/Nx=1/Ny=1/128. We denote ϕxi,yj,tn=ϕi−0.5h,j−0.5h,nΔt by ϕijn for i=1,…,Nx and j=1,…,Ny. An initial condition is given by a simple extension of equation (4):(12)ϕx,y,0=0.1cos2πxonΩ.We can minimize the influence of domain in this way. Therefore, the two-dimensional reference numerical solution is extended from that of the one-dimensional as follows:(13)ϕij=ϕi,for1≤i≤Nx,1≤j≤Ny.Figures5(a) and 5(b) show the initial condition and the reference solution at t=107h4, respectively.Figure 5
(a, b) The initial condition and the reference solution att=107h4, respectively.
(a)(b)We define the discretel2-norm and the discrete maximum norm as(14)e2=1NxNy∑i=1Nx∑j=1Nyeij2,e∞=max1≤i≤Nx1≤j≤Nyeij,respectively. Tables 5 and 6 list the errors and convergence rates of the numerical results at t=107h4. From the results, we observe that both the schemes achieve the expected convergence rates (the first-order for LSS (10) and the second-order for CN scheme (11)).Table 5
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt2500h4Rate1250h4Rate625h4Rate312.5h4e20.00110.9900.00060.9950.00030.9970.0001e∞0.00180.9900.00090.9950.00050.9980.0002Table 6
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 91.9991.600e − 9e∞1.608e − 72.0004.020e − 81.9991.010e − 81.9972.500e − 9


#### 3.2.3. Three-Dimensional Space

Finally, the two numerical schemes are extended to the three-dimensional CH equation as follows:(15)ϕijkn+1−ϕijknΔt=Δhϕijkn3−3ϕijkn+2ϕijkn+1−ε2Δhϕijkn+1,(16)ϕijkn+1−ϕijknΔt=12Δhϕijkn3−ϕijkn−ε2Δhϕijkn+ϕijkn+13−ϕijkn+1−ε2Δhϕijkn+1,where Δhϕijk=ϕi−1,jk+ϕi+1,jk+ϕi,j−1,k+ϕi,j+1,k+ϕij,k−1+ϕij,k+1−6ϕijk/h2. Here, h=1/Nx=1/Ny=1/Nz=1/128. An initial condition is given by(17)ϕx,y,z,0=0.1cos2πxonΩ=0,13.We extend the one-dimensional reference solution to the three-dimensional reference numerical solution:(18)ϕijk=ϕi,for1≤i≤Nx,1≤j≤Ny,1≤k≤Nz.Note that the final timeT is set to 106h4 for ease of computation in the three-dimensional space. Figures 6(a) and 6(b) show the initial condition and the reference solution at t=106h4, respectively.Figure 6
(a, b) The initial condition and the reference solution att=106h4, respectively.
(a)(b)The discretel2-norm and the discrete maximum norm are simply extended as follows:(19)e2=1NxNyNz∑i=1Nx∑j=1Ny∑k=1Nzeijk2,e∞=max1≤i≤Nx1≤j≤Ny1≤k≤Nzeijk,respectively. Tables 7 and 8 list errors and convergence rates for the numerical results at t=106h4. Tables 7 and 8 show that equations (15) and (16) are of first-order and of second-order in time, respectively.Table 7
Errors and convergence rates for the numerical results of the LSS at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e24.287e − 50.9942.153e − 50.9971.079e − 50.9985.400e − 6e∞5.522e − 50.9902.782e − 50.9951.397e − 50.9977.000e − 6Table 8
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.911e − 92.0004.780e − 102.0001.190e − 101.9993.000e − 11e∞2.860e − 92.0007.150e − 101.9981.790e − 101.9554.600e − 11


## 3.1. Convergence and Stability Tests

First, we conduct a traditional convergence test. As a concrete example, an initial condition is given by(4)ϕx,0=0.1cos2πxonx∈Ω=0,1.Note that equation (4) can be used in testing numerical schemes with both the Neumann and periodic boundary conditions. We employ ε=16h/22tanh−10.9≈0.03, which is the approximate transition layer width [27]. Table 1 shows the rate of convergence. We measure the accuracy of time at the final time T using eΔt,0.5Δt2 which is defined as follows:(5)eΔt,Δt/22=1Nx∑i=1NxeiΔt,Δt/22,where eiΔt,Δt/2 is an error between the computed solutions using large and small time steps, respectively. More precisely, it is(6)eiΔt,Δt/2=ϕiΔt−ϕiΔt/2.Table 1
Rate of convergence of RK4.
Δt,Δt/2h4,h4/2Rateh4/2,h4/4Rateh4/4,h4/8eΔt,Δt/224.378e–114.0302.681e–124.0971.567e–13In this test, we adopt the space step size ash=1/10 and the final time as T=10h4. We can confirm that RK4 is of fourth-order in time indeed according to Table 1.For the next test, we investigate the maximum time step sizesΔtc with respect to the spatial step size h. Here, we fix Nt=105 and other parameters and the initial condition are the same as those described above. Table 2 lists the largest time step sizes satisfying that the numerical solution does not blow up.Table 2
Practical maximum time step sizeΔtc with respect to h.
h1/161/321/641/128Δtc4.686e−51.152e−52.810e−67.085e−7Figure1 depicts the practically stable region and maximum time step size Δtc as described in Table 2. The explicit scheme is theoretically Δt≈Oh4.Figure 1
Practical maximum time step sizeΔtc with respect to h.Figure2 represents the CPU times (seconds) with respect to the number of spatial step Nx, performed on Intel Core i5-6400 CPU at 2.70 GHz with 4 GB RAM. When we vary Nx as 16, 32, 64, and 128 with fixed parameters, Δt=1/1284 and Nt=107, the CPU times are 88.812, 175.468, 354.14, and 708.718. That is, when Nx is doubled, the CPU time is approximately doubled.Figure 2
CPU times versus spatial step.


## 3.2. Benchmark Problem

Let us consider the benchmark problem for the CH equation. We adopt equation (4) as an initial condition. For the reference numerical solution of the benchmark problem, we use h=1/128, Δtref=h4, and Nt=107. Therefore, the final time is T=NtΔtref=107h4. For the model parameter, we use ε=24h/22tanh−10.9≈0.045. Figure 3 shows the initial condition and the reference numerical solution at t=107h4.Figure 3
Initial condition and the reference numerical solution att=107h4.


### 3.2.1. One-Dimensional Space

Let us define an errorex,T as the difference between ϕΔt and ϕΔtref at the final time T, i.e., ex,T=ϕΔtx,T−ϕΔtrefx,T, where ϕΔt and ϕΔtref are the numerical solution with a time step Δt and the reference numerical solution with a time step Δtref, respectively. The one-dimensional LSS and CN schemes to be benchmarked are as follows:(7)ϕin+1−ϕinΔt=Δhϕin3−3ϕin+2ϕin+1−ε2Δhϕin+1,(8)ϕin+1−ϕinΔt=12Δhϕin3−ϕin−ε2Δhϕin+ϕin+13−ϕin+1−ε2Δhϕin+1,where Δt is a time step. We solve equations (7) and (8) using a multigrid method [25, 28]. Note that we use the Gauss–Seidel iterative method in multigrid process combined with Newton’s approximation to compute the nonlinear term in equation (8). Figures 4(a) and 4(b) show the errors for the LSS and CN scheme, respectively.Figure 4
Errors between the reference solution and the solutions evaluated using (a) the LSS and (b) the CN scheme att=107h4.
(a)(b)Now, we define the discretel2-norm and the discrete maximum norm as(9)e2=1Nx∑i=1Nxei2,e∞=max1≤i≤Nxei,respectively. Tables 3 and 4 list the numerical errors and convergence rates of the LSS and those of the CN scheme at t=107h4, respectively. We confirm that both the schemes achieve the expected convergence rates (the first-order for the LSS and the second-order for the CN scheme).Table 3
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt1250h4Rate625h4Rate312.5h4Rate156.25h4e25.687e − 40.9952.854e − 40.9971.430e − 40.9997.150e − 5e∞8.970e − 40.9954.500e − 40.9982.254e − 40.9991.128e − 4Table 4
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 92.0011.600e − 9e∞1.608e − 72.0004.020e − 82.0001.010e − 82.0002.500e − 9


### 3.2.2. Two-Dimensional Space

Next, we consider the two-dimensional version of equations (7) and (8) to the CH equation on Ω=0,1×0,1. Straightforward extensions of the LSS and CN schemes are as follows:(10)ϕijn+1−ϕijnΔt=Δhϕijn3−3ϕijn+2ϕijn+1−ε2Δhϕijn+1,(11)ϕijn+1−ϕijnΔt=12Δhϕijn3−ϕijn−ε2Δhϕijn+ϕijn+13−ϕijn+1−ε2Δhϕijn+1,where Δhψij=ψi−1,j+ψi+1,j+ψi,j−1+ψi,j+1−4ψij/h2 and h=1/Nx=1/Ny=1/128. We denote ϕxi,yj,tn=ϕi−0.5h,j−0.5h,nΔt by ϕijn for i=1,…,Nx and j=1,…,Ny. An initial condition is given by a simple extension of equation (4):(12)ϕx,y,0=0.1cos2πxonΩ.We can minimize the influence of domain in this way. Therefore, the two-dimensional reference numerical solution is extended from that of the one-dimensional as follows:(13)ϕij=ϕi,for1≤i≤Nx,1≤j≤Ny.Figures5(a) and 5(b) show the initial condition and the reference solution at t=107h4, respectively.Figure 5
(a, b) The initial condition and the reference solution att=107h4, respectively.
(a)(b)We define the discretel2-norm and the discrete maximum norm as(14)e2=1NxNy∑i=1Nx∑j=1Nyeij2,e∞=max1≤i≤Nx1≤j≤Nyeij,respectively. Tables 5 and 6 list the errors and convergence rates of the numerical results at t=107h4. From the results, we observe that both the schemes achieve the expected convergence rates (the first-order for LSS (10) and the second-order for CN scheme (11)).Table 5
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt2500h4Rate1250h4Rate625h4Rate312.5h4e20.00110.9900.00060.9950.00030.9970.0001e∞0.00180.9900.00090.9950.00050.9980.0002Table 6
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 91.9991.600e − 9e∞1.608e − 72.0004.020e − 81.9991.010e − 81.9972.500e − 9


### 3.2.3. Three-Dimensional Space

Finally, the two numerical schemes are extended to the three-dimensional CH equation as follows:(15)ϕijkn+1−ϕijknΔt=Δhϕijkn3−3ϕijkn+2ϕijkn+1−ε2Δhϕijkn+1,(16)ϕijkn+1−ϕijknΔt=12Δhϕijkn3−ϕijkn−ε2Δhϕijkn+ϕijkn+13−ϕijkn+1−ε2Δhϕijkn+1,where Δhϕijk=ϕi−1,jk+ϕi+1,jk+ϕi,j−1,k+ϕi,j+1,k+ϕij,k−1+ϕij,k+1−6ϕijk/h2. Here, h=1/Nx=1/Ny=1/Nz=1/128. An initial condition is given by(17)ϕx,y,z,0=0.1cos2πxonΩ=0,13.We extend the one-dimensional reference solution to the three-dimensional reference numerical solution:(18)ϕijk=ϕi,for1≤i≤Nx,1≤j≤Ny,1≤k≤Nz.Note that the final timeT is set to 106h4 for ease of computation in the three-dimensional space. Figures 6(a) and 6(b) show the initial condition and the reference solution at t=106h4, respectively.Figure 6
(a, b) The initial condition and the reference solution att=106h4, respectively.
(a)(b)The discretel2-norm and the discrete maximum norm are simply extended as follows:(19)e2=1NxNyNz∑i=1Nx∑j=1Ny∑k=1Nzeijk2,e∞=max1≤i≤Nx1≤j≤Ny1≤k≤Nzeijk,respectively. Tables 7 and 8 list errors and convergence rates for the numerical results at t=106h4. Tables 7 and 8 show that equations (15) and (16) are of first-order and of second-order in time, respectively.Table 7
Errors and convergence rates for the numerical results of the LSS at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e24.287e − 50.9942.153e − 50.9971.079e − 50.9985.400e − 6e∞5.522e − 50.9902.782e − 50.9951.397e − 50.9977.000e − 6Table 8
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.911e − 92.0004.780e − 102.0001.190e − 101.9993.000e − 11e∞2.860e − 92.0007.150e − 101.9981.790e − 101.9554.600e − 11


## 3.2.1. One-Dimensional Space

Let us define an errorex,T as the difference between ϕΔt and ϕΔtref at the final time T, i.e., ex,T=ϕΔtx,T−ϕΔtrefx,T, where ϕΔt and ϕΔtref are the numerical solution with a time step Δt and the reference numerical solution with a time step Δtref, respectively. The one-dimensional LSS and CN schemes to be benchmarked are as follows:(7)ϕin+1−ϕinΔt=Δhϕin3−3ϕin+2ϕin+1−ε2Δhϕin+1,(8)ϕin+1−ϕinΔt=12Δhϕin3−ϕin−ε2Δhϕin+ϕin+13−ϕin+1−ε2Δhϕin+1,where Δt is a time step. We solve equations (7) and (8) using a multigrid method [25, 28]. Note that we use the Gauss–Seidel iterative method in multigrid process combined with Newton’s approximation to compute the nonlinear term in equation (8). Figures 4(a) and 4(b) show the errors for the LSS and CN scheme, respectively.Figure 4
Errors between the reference solution and the solutions evaluated using (a) the LSS and (b) the CN scheme att=107h4.
(a)(b)Now, we define the discretel2-norm and the discrete maximum norm as(9)e2=1Nx∑i=1Nxei2,e∞=max1≤i≤Nxei,respectively. Tables 3 and 4 list the numerical errors and convergence rates of the LSS and those of the CN scheme at t=107h4, respectively. We confirm that both the schemes achieve the expected convergence rates (the first-order for the LSS and the second-order for the CN scheme).Table 3
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt1250h4Rate625h4Rate312.5h4Rate156.25h4e25.687e − 40.9952.854e − 40.9971.430e − 40.9997.150e − 5e∞8.970e − 40.9954.500e − 40.9982.254e − 40.9991.128e − 4Table 4
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 92.0011.600e − 9e∞1.608e − 72.0004.020e − 82.0001.010e − 82.0002.500e − 9


## 3.2.2. Two-Dimensional Space

Next, we consider the two-dimensional version of equations (7) and (8) to the CH equation on Ω=0,1×0,1. Straightforward extensions of the LSS and CN schemes are as follows:(10)ϕijn+1−ϕijnΔt=Δhϕijn3−3ϕijn+2ϕijn+1−ε2Δhϕijn+1,(11)ϕijn+1−ϕijnΔt=12Δhϕijn3−ϕijn−ε2Δhϕijn+ϕijn+13−ϕijn+1−ε2Δhϕijn+1,where Δhψij=ψi−1,j+ψi+1,j+ψi,j−1+ψi,j+1−4ψij/h2 and h=1/Nx=1/Ny=1/128. We denote ϕxi,yj,tn=ϕi−0.5h,j−0.5h,nΔt by ϕijn for i=1,…,Nx and j=1,…,Ny. An initial condition is given by a simple extension of equation (4):(12)ϕx,y,0=0.1cos2πxonΩ.We can minimize the influence of domain in this way. Therefore, the two-dimensional reference numerical solution is extended from that of the one-dimensional as follows:(13)ϕij=ϕi,for1≤i≤Nx,1≤j≤Ny.Figures5(a) and 5(b) show the initial condition and the reference solution at t=107h4, respectively.Figure 5
(a, b) The initial condition and the reference solution att=107h4, respectively.
(a)(b)We define the discretel2-norm and the discrete maximum norm as(14)e2=1NxNy∑i=1Nx∑j=1Nyeij2,e∞=max1≤i≤Nx1≤j≤Nyeij,respectively. Tables 5 and 6 list the errors and convergence rates of the numerical results at t=107h4. From the results, we observe that both the schemes achieve the expected convergence rates (the first-order for LSS (10) and the second-order for CN scheme (11)).Table 5
Errors and convergence rates for the numerical results of the LSS at the final timeT=107h4. Here, h=1/128 is used.
Δt2500h4Rate1250h4Rate625h4Rate312.5h4e20.00110.9900.00060.9950.00030.9970.0001e∞0.00180.9900.00090.9950.00050.9980.0002Table 6
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=107h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.014e − 72.0002.540e − 82.0006.300e − 91.9991.600e − 9e∞1.608e − 72.0004.020e − 81.9991.010e − 81.9972.500e − 9


## 3.2.3. Three-Dimensional Space

Finally, the two numerical schemes are extended to the three-dimensional CH equation as follows:(15)ϕijkn+1−ϕijknΔt=Δhϕijkn3−3ϕijkn+2ϕijkn+1−ε2Δhϕijkn+1,(16)ϕijkn+1−ϕijknΔt=12Δhϕijkn3−ϕijkn−ε2Δhϕijkn+ϕijkn+13−ϕijkn+1−ε2Δhϕijkn+1,where Δhϕijk=ϕi−1,jk+ϕi+1,jk+ϕi,j−1,k+ϕi,j+1,k+ϕij,k−1+ϕij,k+1−6ϕijk/h2. Here, h=1/Nx=1/Ny=1/Nz=1/128. An initial condition is given by(17)ϕx,y,z,0=0.1cos2πxonΩ=0,13.We extend the one-dimensional reference solution to the three-dimensional reference numerical solution:(18)ϕijk=ϕi,for1≤i≤Nx,1≤j≤Ny,1≤k≤Nz.Note that the final timeT is set to 106h4 for ease of computation in the three-dimensional space. Figures 6(a) and 6(b) show the initial condition and the reference solution at t=106h4, respectively.Figure 6
(a, b) The initial condition and the reference solution att=106h4, respectively.
(a)(b)The discretel2-norm and the discrete maximum norm are simply extended as follows:(19)e2=1NxNyNz∑i=1Nx∑j=1Ny∑k=1Nzeijk2,e∞=max1≤i≤Nx1≤j≤Ny1≤k≤Nzeijk,respectively. Tables 7 and 8 list errors and convergence rates for the numerical results at t=106h4. Tables 7 and 8 show that equations (15) and (16) are of first-order and of second-order in time, respectively.Table 7
Errors and convergence rates for the numerical results of the LSS at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e24.287e − 50.9942.153e − 50.9971.079e − 50.9985.400e − 6e∞5.522e − 50.9902.782e − 50.9951.397e − 50.9977.000e − 6Table 8
Errors and convergence rates for the numerical results of the CN scheme at the final timeT=106h4. Here, h=1/128 is used.
Δt10000h4Rate5000h4Rate2500h4Rate1250h4e21.911e − 92.0004.780e − 102.0001.190e − 101.9993.000e − 11e∞2.860e − 92.0007.150e − 101.9981.790e − 101.9554.600e − 11


## 4. Conclusions

In this work, we proposed a very simple benchmark problem, which is a cosine function as the initial condition for the numerical methods of the CH equations. For the numerical solution, we used RK4 for the temporal integration. Using the proposed benchmark problem solution, we performed the convergence tests for LSS and CN schemes and obtained the expected first-order and second-order convergence rates in the one-, two-, and three-dimensional CH equations, respectively. The proposed approach is less biased and provides more rational results compared to conventional self-convergence tests. In future work, we will design simple benchmark problems for other nonlinear partial differential equations such as the nonlocal CH, AC, and nonlocal AC equations.



---

*Source: 8889603-2021-03-09.xml*