Next Article in Journal
On Formality of Some Homogeneous Spaces
Previous Article in Journal
Cosmological Consequences of a Parametrized Equation of State
Previous Article in Special Issue
A Numerical Solution of Fredholm Integral Equations of the Second Kind Based on Tight Framelets Generated by the Oblique Extension Principle
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Simulation of Pattern Formation on Surfaces Using an Efficient Linear Second-Order Method

Department of Mathematics, Kwangwoon University, Seoul 01897, Korea
Symmetry 2019, 11(8), 1010; https://doi.org/10.3390/sym11081010
Submission received: 4 June 2019 / Revised: 24 July 2019 / Accepted: 25 July 2019 / Published: 5 August 2019
(This article belongs to the Special Issue Numerical Analysis or Numerical Method in Symmetry)

Abstract

:
We present an efficient linear second-order method for a Swift–Hohenberg (SH) type of a partial differential equation having quadratic-cubic nonlinearity on surfaces to simulate pattern formation on surfaces numerically. The equation is symmetric under a change of sign of the density field if there is no quadratic nonlinearity. We introduce a narrow band neighborhood of a surface and extend the equation on the surface to the narrow band domain. By applying a pseudo-Neumann boundary condition through the closest point, the Laplace–Beltrami operator can be replaced by the standard Laplacian operator. The equation on the narrow band domain is split into one linear and two nonlinear subequations, where the nonlinear subequations are independent of spatial derivatives and thus are ordinary differential equations and have closed-form solutions. Therefore, we only solve the linear subequation on the narrow band domain using the Crank–Nicolson method. Numerical experiments on various surfaces are given verifying the accuracy and efficiency of the proposed method.

1. Introduction

A Swift–Hohenberg (SH) type of partial differential equation [1] has been used to study pattern formation [2,3,4,5]:
ϕ t = ϕ 3 g ϕ 2 + ϵ + ( 1 + Δ ) 2 ϕ ,
where ϕ is the density field and g 0 and ϵ > 0 are constants. In general, the equation does not have an analytical solution, thus various computational algorithms [6,7,8,9,10,11,12,13] have been proposed to obtain a numerical solution. However, most of them were solved on flat surfaces except [12,13].
In this paper, we present an efficient linear second-order method for the SH type of equation on surfaces, which is based on the closest point method [14,15]. We introduce a narrow band domain of a surface and apply a pseudo-Neumann boundary condition on the boundary of the narrow band domain through the closest point [16]. This results in a constant value of ϕ in the direction normal to the surface, thus the Laplace–Beltrami operator can be replaced by the standard Laplacian operator. In addition, we split the equation into one linear and two nonlinear subequations [17,18], where the nonlinear subequations are independent of spatial derivatives and thus are ordinary differential equations and have closed-form solutions. Therefore, we only solve the linear subequation on the narrow band domain using the Crank–Nicolson method. As a result, our method is easy to implement and linear.
This paper is organized as follows. In Section 2, we describe the SH type of equation on a narrow band domain. In Section 3, we propose an efficient linear second-order method for the equation on the narrow band domain. Numerical examples on various surfaces are given in Section 4. Finally, we conclude in Section 5.

2. Swift–Hohenberg Type of Equation on a Narrow Band Domain

The SH type of equation on a surface S is given by
ϕ ( x , t ) t = ϕ 3 ( x , t ) g ϕ 2 ( x , t ) + ϵ + ( 1 + Δ S ) 2 ϕ ( x , t ) , x S , 0 < t T ,
where Δ S is the Laplace–Beltrami operator [19,20]. Next, let Ω δ = { y | x S , y = x + η n ( x ) for | η | < δ } be a δ -neighborhood of S , where n ( x ) is a unit normal vector at x . Then, we extend the Equation (1) to the narrow band domain Ω δ :
ϕ ( x , t ) t = ϕ 3 ( x , t ) g ϕ 2 ( x , t ) + ϵ + ( 1 + Δ S ) 2 ϕ ( x , t ) , x Ω δ , 0 < t T
with the pseudo-Neumann boundary condition on Ω δ :
ϕ ( x , t ) = ϕ ( cp ( x ) , t ) ,
where cp ( x ) is a point on S , which is closest to x Ω δ [14]. For a sufficiently small δ , ϕ is constant in the direction normal to the surface. Thus, the Laplace–Beltrami operator in Ω δ can be replaced by the standard Laplacian operator [14], i.e.,
ϕ ( x , t ) t = ϕ 3 ( x , t ) g ϕ 2 ( x , t ) + ϵ + ( 1 + Δ ) 2 ϕ ( x , t ) , x Ω δ , 0 < t T .

3. Numerical Method

In this section, we propose an efficient linear second-order method for solving Equation (4) with the boundary condition (3). We discretize Equation (4) in Ω = [ L x / 2 , L x / 2 ] × [ L y / 2 , L y / 2 ] × [ L z / 2 , L z / 2 ] that includes Ω δ . Let h = L x / N x = L y / N y = L z / N z be the uniform grid size, where N x , N y , and N z are positive integers. Let Ω h = { x i j k = ( x i , y j , z k ) | x i = L x / 2 + i h , y j = L y / 2 + j h , z k = L z / 2 + k h for 0 i N x , 0 j N y , 0 k N z } be a discrete domain. Let ϕ i j k n be an approximation of ϕ ( x i j k , n Δ t ) , where Δ t is the time step. Let Ω δ h = { x i j k | | ψ i j k | < δ } be a discrete narrow band domain, where ψ is a signed distance function for the surface S , and Ω δ h = { x i j k | I i j k | h I i j k | 0 } are discrete domain boundary points, where h I i j k = ( I i + 1 , j , k I i 1 , j , k , I i , j + 1 , k I i , j 1 , k , I i , j , k + 1 I i , j , k 1 ) / ( 2 h ) . Here, I i j k = 0 if x i j k Ω δ h , and I i j k = 1 , otherwise.
We here split Equation (4) into the following subequations:
ϕ t = ( ϕ 3 ϵ ϕ ) ,
ϕ t = g ϕ 2 ,
ϕ t = ( 1 + Δ ) 2 ϕ .
Equations (5) and (6) are solved analytically and the solutions ϕ i j k n + 1 are given as follows:
ϕ i j k n + 1 = ϕ i j k n ( ϕ i j k n ) 2 / ϵ + ( 1 ( ϕ i j k n ) 2 / ϵ ) e 2 ϵ Δ t and ϕ i j k n + 1 = ϕ i j k n 1 g Δ t ϕ i j k n ,
respectively. In addition, Equation (7) is solved using the Crank–Nicolson method:
ϕ i j k n + 1 ϕ i j k n Δ t = ( 1 + Δ h ) 2 2 ( ϕ i j k n + 1 + ϕ i j k n )
with the boundary condition on Ω δ h :
ϕ i j k n = ϕ n ( cp ( x i j k ) ) .
Here, Δ h ϕ i j k = ( ϕ i + 1 , j , k + ϕ i 1 , j , k + ϕ i , j + 1 , k + ϕ i , j 1 , k + ϕ i , j , k + 1 + ϕ i , j , k 1 6 ϕ i j k ) / h 2 . The numerical closest point cp ( x i j k ) for a point x i j k Ω δ h is defined as
cp ( x i j k ) = x i j k | ψ i j k | h | ψ i j k | | h | ψ i j k | | .
In general, cp ( x i j k ) is not a grid point in Ω δ h , i.e., cp ( x i j k ) Ω δ h , and thus we use trilinear interpolation and take δ > 3 h to obtain ϕ ( cp ( x i j k ) ) . The resulting implicit linear discrete system of Equation (8) is solved efficiently using the Jacobi iterative method. We iterate the Jacobi iteration until a discrete L 2 -norm of the consecutive error on Ω δ h is less than a tolerance t o l . Here, the discrete L 2 -norm on Ω δ h is defined as ϕ L 2 ( Ω δ h ) = x i j k Ω δ h ϕ i j k 2 / # Ω δ h , where # Ω δ h is the cardinality of Ω δ h . Then, the second-order solution of Equation (4) is evolved by five stages [21]
ϕ i j k ( 1 ) = ϕ i j k n ( ϕ i j k n ) 2 / ϵ + ( 1 ( ϕ i j k n ) 2 / ϵ ) e ϵ Δ t , ϕ i j k ( 2 ) = ϕ i j k ( 1 ) 1 ( g Δ t / 2 ) ϕ i j k ( 1 ) , ϕ i j k ( 3 ) ϕ i j k ( 2 ) Δ t = ( 1 + Δ h ) 2 2 ( ϕ i j k ( 3 ) + ϕ i j k ( 2 ) ) , ϕ i j k ( 4 ) = ϕ i j k ( 3 ) 1 ( g Δ t / 2 ) ϕ i j k ( 3 ) , ϕ i j k n + 1 = ϕ i j k ( 4 ) ( ϕ i j k ( 4 ) ) 2 / ϵ + ( 1 ( ϕ i j k ( 4 ) ) 2 / ϵ ) e ϵ Δ t .

4. Numerical Experiments

4.1. Convergence Test

In order to verify the rate of convergence of the proposed method, we consider the evolution of ϕ on a unit sphere. An initial piece of data is
ϕ ( x , y , z , 0 ) = 0.15 + 0.1 cos ( 2 π x ) cos ( 2 π y ) cos ( 2 π z )
and a signed distance function for the unit sphere is
ψ ( x , y , z ) = x 2 + y 2 + z 2 1
on Ω = [ 1 . 5 , 1 . 5 ] 3 . We fix the grid size to h = 0.125 and vary Δ t = T / 2 , T / 2 2 , T / 2 3 , T / 2 4 for T = 0.00025 with ϵ = 0.25 , δ = 2.2 3 h , and t o l = Δ t . Table 1 shows the L 2 -errors of ϕ ( x , y , z , T ) and convergence rates with g = 0 . Here, the errors are computed by comparison with a reference numerical solution using Δ t = T / 2 6 . It is observed that the method is second-order accurate in time. Note that we obtain the same result for g = 1 .

4.2. Pattern Formation on a Sphere

Unless otherwise stated, we take an initial piece of data as
ϕ ( x , y , z , 0 ) = 0.15 + rand ( x , y , z ) ,
where rand ( x , y , z ) is a uniformly distributed random number between 0.1 and 0.1 at the grid points, and use ϵ = 0.25 , h = 1 , Δ t = 0.1 , δ = 1.1 3 h , and t o l = 10 4 .
For g = 0 and 1, Figure 1 and Figure 2 show the evolution of ϕ ( x , y , z , t ) on a sphere with ψ ( x , y , z ) = x 2 + y 2 + z 2 32 on Ω = [ 36 , 36 ] 3 , respectively. Depending on the value of g, we have different patterns, such as striped (Figure 1) and hexagonal (Figure 2) [11]. Figure 3 shows the energy decay with g = 0 and 1, where the energy E ( ϕ ) is defined by
E ( ϕ ) = Ω δ 1 4 ϕ 4 g 3 ϕ 3 + 1 2 ϕ ϵ + ( 1 + Δ ) 2 ϕ d x .

4.3. Pattern Formation on a Sphere Perturbed by a Spherical Harmonic

In this section, we perform the evolution of ϕ on a sphere of center ( 0 , 0 , 0 ) and radius 32 perturbed by a spherical harmonic 10 Y 10 7 ( θ , φ ) . Here, θ and φ are the polar and azimuthal angles, respectively, and the computational domain is Ω = [ 40 , 40 ] 3 . Figure 4 and Figure 5 show the evolution of ϕ ( x , y , z , t ) with g = 0 and 1, respectively. From the results in Figure 4 and Figure 5, we can see that our method can solve the SH type of equation on not only simple but also complex surfaces. Figure 6 shows the energy decay with g = 0 and 1.

4.4. Pattern Formation on a Spindle

Finally, we simulate the evolution of ϕ on a spindle that has narrow and sharp tips. The spindle is defined parametrically as
x = 16 cos θ sin φ , y = 16 sin θ sin φ , z = 32 2 φ π 1 ,
where θ [ 0 , 2 π ) and φ [ 0 , π ) , and the computational domain is Ω = [ 20 , 20 ] × [ 20 , 20 ] × [ 36 , 36 ] . Figure 7 and Figure 8 show the evolution of ϕ ( x , y , z , t ) with g = 0 and 1, respectively. The results in Figure 7 and Figure 8 suggest that pattern formation on a surface having narrow and sharp tips can be simulated by using our method. Figure 9 shows the energy decay with g = 0 and 1.

5. Conclusions

We simulated pattern formation on surfaces numerically by solving the SH type of equation on surfaces by using the efficient linear second-order method. The method was based on the closest point and operator splitting methods and thus was easy to implement and linear. We confirmed that the proposed method gives the desired order of accuracy in time and observed that pattern formation on surfaces is affected by the value of g.

Funding

This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2019R1C1C1011112).

Acknowledgments

The corresponding author thanks the reviewers for the constructive and helpful comments on the revision of this article.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Swift, J.; Hohenberg, P.C. Hydrodynamic fluctuations at the convective instability. Phys. Rev. A 1977, 15, 319–328. [Google Scholar] [CrossRef] [Green Version]
  2. Hohenberg, P.C.; Swift, J.B. Effects of additive noise at the onset of Rayleigh–Bénard convection. Phys. Rev. A 1992, 46, 4773–4785. [Google Scholar] [CrossRef] [PubMed]
  3. Cross, M.C.; Hohenberg, P.C. Pattern formation outside of equilibrium. Rev. Mod. Phys. 1993, 65, 851–1112. [Google Scholar] [CrossRef] [Green Version]
  4. Rosa, R.R.; Pontes, J.; Christov, C.I.; Ramos, F.M.; Rodrigues Neto, C.; Rempel, E.L.; Walgraef, D. Gradient pattern analysis of Swift–Hohenberg dynamics: Phase disorder characterization. Phys. A 2000, 283, 156–159. [Google Scholar] [CrossRef]
  5. Hutt, A.; Atay, F.M. Analysis of nonlocal neural fields for both general and gamma-distributed connectivities. Phys. D 2005, 203, 30–54. [Google Scholar] [CrossRef]
  6. Cheng, M.; Warren, J.A. An efficient algorithm for solving the phase field crystal model. J. Comput. Phys. 2008, 227, 6241–6248. [Google Scholar] [CrossRef] [Green Version]
  7. Wise, S.M.; Wang, C.; Lowengrub, J.S. An energy-stable and convergent finite-difference scheme for the phase field crystal equation. SIAM J. Numer. Anal. 2009, 47, 2269–2288. [Google Scholar] [CrossRef]
  8. Gomez, H.; Nogueira, X. A new space–time discretization for the Swift–Hohenberg equation that strictly respects the Lyapunov functional. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 4930–4946. [Google Scholar] [CrossRef]
  9. Elsey, M.; Wirth, B. A simple and efficient scheme for phase field crystal simulation. ESAIM: Math. Model. Numer. Anal. 2013, 47, 1413–1432. [Google Scholar] [CrossRef] [Green Version]
  10. Sarmiento, A.F.; Espath, L.F.R.; Vignal, P.; Dalcin, L.; Parsani, M.; Calo, V.M. An energy-stable generalized-α method for the Swift–Hohenberg equation. J. Comput. Appl. Math. 2018, 344, 836–851. [Google Scholar] [CrossRef]
  11. Lee, H.G. An energy stable method for the Swift–Hohenberg equation with quadratic–cubic nonlinearity. Comput. Methods Appl. Mech. Eng. 2019, 343, 40–51. [Google Scholar] [CrossRef]
  12. Matthews, P.C. Pattern formation on a sphere. Phys. Rev. E 2003, 67, 036206. [Google Scholar] [CrossRef] [PubMed]
  13. Sigrist, R.; Matthews, P. Symmetric spiral patterns on spheres. SIAM J. Appl. Dyn. Syst. 2011, 10, 1177–1211. [Google Scholar] [CrossRef]
  14. Ruuth, S.J.; Merriman, B. A simple embedding method for solving partial differential equations on surfaces. J. Comput. Phys. 2008, 227, 1943–1961. [Google Scholar] [CrossRef]
  15. Choi, Y.; Jeong, D.; Lee, S.; Yoo, M.; Kim, J. Motion by mean curvature of curves on surfaces using the Allen–Cahn equation. Int. J. Eng. Sci. 2015, 97, 126–132. [Google Scholar] [CrossRef]
  16. Lee, H.G.; Kim, J. A simple and efficient finite difference method for the phase-field crystal equation on curved surfaces. Comput. Methods Appl. Mech. Eng. 2016, 307, 32–43. [Google Scholar] [CrossRef]
  17. Pak, D.; Han, C.; Hong, W.-T. Iterative speedup by utilizing symmetric data in pricing options with two risky assets. Symmetry 2017, 9, 12. [Google Scholar] [CrossRef]
  18. Zong, C.; Tang, Y.; Cho, Y.J. Convergence analysis of an inexact three-operator splitting algorithm. Symmetry 2018, 10, 563. [Google Scholar] [CrossRef]
  19. Bertalmío, M.; Cheng, L.-T.; Osher, S.; Sapiro, G. Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. 2001, 174, 759–780. [Google Scholar] [CrossRef]
  20. Greer, J.B. An improvement of a recent Eulerian method for solving PDEs on general geometries. J. Sci. Comput. 2006, 29, 321–352. [Google Scholar] [CrossRef]
  21. Strang, G. On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 1968, 5, 506–517. [Google Scholar] [CrossRef]
Figure 1. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.7540 and 0.7783 , respectively.
Figure 1. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.7540 and 0.7783 , respectively.
Symmetry 11 01010 g001
Figure 2. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.4320 and 0.7152 , respectively.
Figure 2. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.4320 and 0.7152 , respectively.
Symmetry 11 01010 g002
Figure 3. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the sphere with g = 0 and 1.
Figure 3. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the sphere with g = 0 and 1.
Symmetry 11 01010 g003
Figure 4. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.8717 and 0.8372 , respectively.
Figure 4. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.8717 and 0.8372 , respectively.
Symmetry 11 01010 g004
Figure 5. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.4833 and 0.7135 , respectively.
Figure 5. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.4833 and 0.7135 , respectively.
Symmetry 11 01010 g005
Figure 6. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the perturbed sphere with g = 0 and 1.
Figure 6. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the perturbed sphere with g = 0 and 1.
Symmetry 11 01010 g006
Figure 7. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.7059 and 0.7593 , respectively.
Figure 7. Evolution of ϕ ( x , y , z , t ) with g = 0 . The yellow and blue regions indicate ϕ = 0.7059 and 0.7593 , respectively.
Symmetry 11 01010 g007
Figure 8. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.3842 and 0.6224 , respectively.
Figure 8. Evolution of ϕ ( x , y , z , t ) with g = 1 . The yellow and blue regions indicate ϕ = 1.3842 and 0.6224 , respectively.
Symmetry 11 01010 g008
Figure 9. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the spindle with g = 0 and 1.
Figure 9. Evolution of E ( ϕ ) / E ( ϕ 0 ) on the spindle with g = 0 and 1.
Symmetry 11 01010 g009
Table 1. L 2 -errors and convergence rates for g = 0 .
Table 1. L 2 -errors and convergence rates for g = 0 .
Δ t T / 2 T / 2 2 T / 2 3 T / 2 4
L 2 -error 5.445 × 10 3 1.341 × 10 3 2.862 × 10 4 5.519 × 10 5
Rate 2.02 2.22 2.37

Share and Cite

MDPI and ACS Style

Lee, H.G. Numerical Simulation of Pattern Formation on Surfaces Using an Efficient Linear Second-Order Method. Symmetry 2019, 11, 1010. https://doi.org/10.3390/sym11081010

AMA Style

Lee HG. Numerical Simulation of Pattern Formation on Surfaces Using an Efficient Linear Second-Order Method. Symmetry. 2019; 11(8):1010. https://doi.org/10.3390/sym11081010

Chicago/Turabian Style

Lee, Hyun Geun. 2019. "Numerical Simulation of Pattern Formation on Surfaces Using an Efficient Linear Second-Order Method" Symmetry 11, no. 8: 1010. https://doi.org/10.3390/sym11081010

APA Style

Lee, H. G. (2019). Numerical Simulation of Pattern Formation on Surfaces Using an Efficient Linear Second-Order Method. Symmetry, 11(8), 1010. https://doi.org/10.3390/sym11081010

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop