Source-linked AI summary
An Energy Stable BDF2 Fourier Pseudo-Spectral Numerical Scheme for the Square Phase Field Crystal Equation
Kelong Cheng, Cheng Wang, Steven M. Wise
TL;DR
SPFC numerical analysis is challenging because its modified energy yields a highly nonlinear operator composition, and rigorous numerical analysis had remained open. The paper develops second-order energy-stable Fourier pseudo-spectral schemes and analyzes their solvability and convergence. The analysis establishes unique solvability, mass conservation, and an optimal-rate convergence result under the stated setting.
Problem
The highly nonlinear 4-Laplacian composition makes SPFC analysis and numerical approximation especially challenging, and rigorous numerical analysis had remained an open problem.
Source: §1 · p. 3Method
The paper proposes second-order energy-stable Fourier pseudo-spectral schemes using a modified BDF approach with explicit extrapolation of the concave diffusion term.
Source: §1 · p. 3, §7 · p. 20Results
The schemes have unique solutions and conserve mass, and the analysis derives an optimal-rate error bound in the ℓ∞(0,T; ℓ2) ∩ ℓ2(0,T; H_N^3) norm.
Source: §3 · p. 10, §7 · p. 20Takeaways & Limitations
Aliasing-error control and a discrete H2-to-W1,6 embedding provide the nonlinear estimate used to derive the optimal-rate error analysis.
Source: §7 · p. 20The analysis assumes periodic fields on a rectangular domain in D=2 or 3.
Source: §1 · p. 2
Abstract
from arXiv · showhide
In this paper we propose and analyze an energy stable numerical scheme for the square phase field crystal (SPFC) equation, a gradient flow modeling crystal dynamics at the atomic scale in space but on diffusive scales in time. In particular, a modification of the free energy potential to the standard phase field crystal model leads to a composition of the 4-Laplacian and the regular Laplacian operators. To overcome the difficulties associated with this highly nonlinear operator, we design numerical algorithms based on the structures of the individual energy terms. A Fourier pseudo-spectral approximation is taken in space, in such a way that the energy structure is respected, and summation-by-parts formulae enable us to study the discrete energy stability for such a high-order spatial discretization. In the temporal approximation, a second order BDF stencil is applied, combined with an appropriate extrapolation for the concave diffusion term(s). A second order artificial Douglas-Dupont-type regularization term is added to ensure energy stability, and a careful analysis leads to the artificial linear diffusion coming at an order lower that that of surface diffusion term. Such a choice leads to reduced numerical dissipation. At a theoretical level, the unique solvability, energy stability are established, and an optimal rate convergence analysis is derived in the $\ell^\infty (0,T; \ell^2) \cap \ell^2 (0,T; H_N^3)$ norm. In the numerical implementation, the preconditioned steepest descent (PSD) iteration is applied to solve for the composition of the highly nonlinear 4-Laplacian term and the standard Laplacian term, and a geometric convergence is assured for such an iteration. Finally, a few numerical experiments are presented, which confirm the robustness and accuracy of the proposed scheme.
1 Introduction
The SPFC model is difficult to analyze numerically because its modified energy produces a highly nonlinear 4-Laplacian composition, and rigorous analysis had remained open. The paper proposes energy-stable, second-order Fourier pseudo-spectral schemes and analyzes their solvability, stability, convergence, and implementation.
1 Introduction
- Rigorous numerical analysis of the SPFC model remained an open problem because finite-difference multigrid performed poorly on the highly nonlinear 4-Laplacian, whose H−1 gradient flow further complicates solving.An efficient spectral/pseudo-spectral PSD solver for the composition of 4-Laplacian and regular Laplacian operators had not been available in the existing literature.Source: §1 · pp. 3–4
- Two energy-stable SPFC schemes combine second-order temporal accuracy with Fourier pseudo-spectral spatial discretization.They differ in whether the concave term −|∇φ|^2 or −εφ^2 is treated as destabilizing.Source: §1 · p. 3
- Explicit extrapolation of the concave diffusion term in a modified BDF approach ensures unique solvability without sacrificing energy stability.A second-order Douglas-Dupont-type regularization is added, with energy stability guaranteed under a mild condition.Source: §1 · pp. 3–4
- Summation-by-parts formulas and aliasing-error control support unique solvability and energy stability for the fully discrete Fourier pseudo-spectral scheme.The discrete stability yields a uniform-in-time H2 bound, used to obtain a W1,6 bound and an optimal-rate error analysis.Source: §1 · p. 4
- The convergence analysis establishes an optimal O(Δt^2 + h^m) rate in the ℓ∞(0,T; ℓ2) ∩ ℓ2(0,T; H_N^3) norm.The argument uses a discrete Sobolev embedding in the Fourier pseudo-spectral space and aliasing-error estimates.Source: §1 · p. 4
- A preconditioned steepest descent solver addresses the composition of the 4-Laplacian and regular Laplacian, with geometric convergence derived for the iteration.Each search direction requires solving a linear elliptic equation, and the optimization along that direction is one-dimensional.Source: §1 · p. 4
2 The numerical scheme
The numerical scheme is built on a Fourier pseudo-spectral grid representation, with discrete operators, norms, and summation-by-parts identities used in its analysis.
2.1 Fourier pseudo-spectral approximations
- The Fourier pseudo-spectral method represents functions on a quadrature grid, simplifying some operator evaluations and enabling fast FFT-based calculations.The analysis assumes a three-dimensional unit cube, an odd number of grid points per direction, and a uniform grid; the description can be modified for two dimensions.Source: §2.1 · p. 5
- Collocation Fourier derivatives are used to define the discrete gradient, divergence, and Laplacian operators.The operators are defined on the numerical grid.Source: §2.1 · pp. 5–6
- The scheme defines spectral interpolation, grid projection, de-aliasing, and standard Fourier projection operators for representing grid functions and periodic functions.The Fourier projection and de-aliasing operators generally differ, while de-aliasing acts as the identity on functions already in the trigonometric-polynomial space.Source: §2.1 · p. 6
- Summation-by-parts identities provide the discrete operator relations used in the analysis, including for the discrete H−1 norm.The discrete H−1 norm is defined on the mean-zero grid-function space because the SPFC equation is an H−1 gradient flow.Source: §2.1 · p. 7
- The discretization introduces grid-function norms and a discrete SPFC energy for use in the subsequent analysis.The definitions include ℓp norms and discrete H1 and H2 norms.Source: §2.1 · p. 8
3 Solvability and stability of the fully discrete schemes
The two fully discrete schemes are uniquely solvable, mass conservative, and energy dissipative under their stated conditions. Their modified BDF2 construction treats concave diffusion explicitly and adds regularization; energy stability also yields uniform-in-time spatial bounds.
3 Solvability and stability of the fully discrete schemes
- The BDF2 schemes explicitly update concave diffusion terms and add second-order regularization to secure unique solvability and energy stability.Scheme (3.3) uses a different convexity view, treating −εφ² as destabilizing.Source: §3 · p. 9
- Both schemes have unique solutions and conserve the mean value at every time step.The solvability result assumes the two preceding time levels have equal mean.Source: §3 · p. 10
- The respective modified discrete energies dissipate for schemes (3.2) and (3.3), subject to the stated stability conditions.The energy-stability analysis uses an added second-order Douglas-Dupont-type regularization.Source: §3 · pp. 9, 10, 12
- Energy stability supplies a uniform-in-time H_N^2 bound, which combines with discrete Sobolev embedding to give a uniform W_N^{1,6} estimate.The stated H_N^2 bound is independent of the mesh size, time step, and final time.Source: §3 · pp. 13–14
4 Convergence analysis
The convergence analysis establishes an error estimate for sufficiently regular exact solutions, with a condition on the regularization parameter. Its consistency analysis gives second-order temporal and spectral spatial accuracy.
4 Convergence analysis
- The proof covers scheme (3.2); the authors state that the proof for scheme (3.3) is similar.Source: §4 · p. 14
- For exact solutions in regularity class R and A ≥ ε², the convergence estimate has a constant independent of ∆t and h.The theorem assumes sufficient regularity of the exact SPFC solution.Source: §4 · p. 14
- The local truncation error is bounded by O(∆t² + h^m), giving second-order temporal and spectral spatial consistency.The truncation estimate is stated for exact solutions in regularity class R.Source: §4 · p. 14
5 Preconditioned steepest descent solver
The solver recasts each time-level nonlinear system as minimization of a strictly convex discrete energy and applies preconditioned steepest descent. A linearized operator sets the search direction, and the iteration has geometric convergence independent of N.
5 Preconditioned steepest descent solver
- At each time level, scheme (3.2) is solved as minimization of a strictly convex discrete energy.The fully discrete scheme is the energy's discrete variation set equal to zero.Source: §5 · p. 17
- The PSD iteration uses a linearized nonlinear operator as a preconditioner to choose the steepest-descent search direction.The search-direction problem can be solved efficiently with the FFT.Source: §5 · p. 17
- The iteration converges geometrically at a rate independent of N.The convergence statement compares iterates with the exact numerical solution at the current time level.Source: §5 · p. 17
6 Numerical results
The numerical tests show second-order temporal accuracy and spatial spectral accuracy, while simulations exhibit square-symmetry patterns and energy dissipation. The two temporal schemes retain second-order accuracy, with BDF2-ES-2 producing smaller errors in the reported test.
6.1 Convergence test for the numerical scheme
- The spatial tests show spectral accuracy for the phase variable, with accuracy saturation at increasing N due to round-off errors.The spatial-error comparison fixes Δt = 10^-4 and absorbs temporal errors into the external source term.Source: §6.1 · pp. 18–19
- Both schemes show second-order temporal accuracy, and BDF2-ES-2 has lower error than BDF2-ES-1 in the reported test.The comparison fixes N = 128 and uses time steps from T/100 to T/800 at T = 0.16.Source: §6.1 · p. 19, Figure 2 · p. 19
6.2 Numerical simulation of square symmetry patterns
- Simulations with magnitude-10 nucleation sites at (50,50) or at (25,25), (25,75), (75,25), and (75,75) confirm emergence of rotationally invariant square-symmetry patterns in the density field.The time step is 0.05 on [0, 1000] and 0.1 on [1000, 9000].Source: §6.2 · p. 20
- For the one-site simulation, the energy dissipates and the results for Δt = 0.1, 0.05, and 0.025 agree closely through t = 1000.The energy plots overlap so their differences are indistinguishable.Source: §6.2 · p. 20
7 Concluding remarks
The concluding remarks describe two energy-stable Fourier pseudo-spectral schemes and their analysis for the SPFC equation. They report an optimal-rate error analysis and a geometrically convergent PSD iteration for the composed nonlinear operator.
7 Concluding remarks
- The schemes use a second-order BDF stencil, extrapolation for the concave diffusion term, and lower-order Douglas-Dupont regularization to establish unique solvability and energy stability.The regularization is lower order than surface diffusion, avoiding higher-order artificial diffusion and leading to expected reduced numerical dissipation.Source: §7 · p. 20
- An aliasing-control technique enables a discrete H2-to-W1,6 embedding, yielding optimal-rate error analysis in ℓ∞(0,T;ℓ2) ∩ ℓ2(0,T;H_N^3).The energy stability also yields a uniform-in-time H2 bound for the numerical solution.Source: §7 · p. 20
- A preconditioned steepest descent iteration for the composed 4-Laplacian and regular Laplacian has proven geometric convergence; simulations demonstrate square patterns and scheme accuracy.The paper identifies this as the first such 4-Laplacian solver result in an H−1 gradient flow.Source: §7 · p. 20
A Proof of Proposition 2.7
The proof represents discrete derivatives through the spectral extension of grid functions, then uses continuous Sobolev embedding to bound the derivative estimate.
A Proof of Proposition 2.7
- Periodic boundary conditions and cell-centered representation give the grid function a discrete Fourier transform, which is extended to a continuous function.The proof focuses on estimating the x-derivative as part of the estimate for ∥∇_N f∥_6.Source: A Proof of Proposition 2.7 · p. 21
- Writing g = D_x f pointwise identifies its spectral representation as g_N = ∂_x f_N.The proof applies a preliminary estimate to the trigonometric polynomial representing f.Source: A Proof of Proposition 2.7 · p. 21
- The derivative bound uses the continuous Sobolev embedding ∥∂_x f_N∥_L6 ≤ C∥Δf_N∥_L2 and the unique correspondence between f and f_N.The same type of estimate is stated to hold for the other coordinate directions.Source: A Proof of Proposition 2.7 · p. 21