Convergence Analysis of a Preconditioned Steepest Descent Solver for the Cahn-Hilliard Equation with Logarithmic Potential

TL;DR

This work tackles the numerical solution of the Cahn–Hilliard equation with the physically realistic Flory–Huggins logarithmic potential, which introduces singular behavior near $\phi=\pm1$. It combines a first-order convex-splitting scheme with a preconditioned steepest descent (PSD) iterative solver, using a linearized preconditioner $A_h$ to stabilize and accelerate convergence. The authors prove that the PSD iterations preserve the positivity of the logarithmic arguments at every step and establish a geometric convergence rate under a strict separation from the singular limits, with the rate governed by constants that depend on $\Delta t$, $\varepsilon$, and the separation $\epsilon_0$. Numerical experiments in 2D demonstrate robust performance: rapid convergence (5–10 iterations per time step) and an energy decay consistent with the expected asymptotics, along with a coarsening dynamics that respects the separation property. Overall, the PSD approach provides a fast, reliable nonlinear solver for a challenging gradient-flow problem with a singular energy potential and opens avenues for extensions to 3D and different boundary conditions.

Abstract

In this paper, we provide a theoretical analysis for a preconditioned steepest descent (PSD) iterative solver that improves the computational time of a finite difference numerical scheme for the Cahn-Hilliard equation with Flory-Huggins energy potential. In the numerical design, a convex splitting approach is applied to the chemical potential such that the logarithmic and the surface diffusion terms are treated implicitly while the expansive concave term is treated with an explicit update. The nonlinear and singular nature of the logarithmic energy potential makes the numerical implementation very challenging. However, the positivity-preserving property for the logarithmic arguments, unconditional energy stability, and optimal rate error estimates have been established in a recent work and it has been shown that successful solvers ensure a similar positivity-preserving property at each iteration stage. Therefore, in this work, we will show that the PSD solver ensures a positivity-preserving property at each iteration stage. The PSD solver consists of first computing a search direction (involved with solving a Poisson-like equation) and then takes a one-parameter optimization step over the search direction in which the Newton iteration becomes very powerful. A theoretical analysis is applied to the PSD iteration solver and a geometric convergence rate is proved for the iteration. In particular, the strict separation property of the numerical solution, which indicates a uniform distance between the numerical solution and the singular limit values of $\pm 1$ for the phase variable, plays an essential role in the iteration convergence analysis. A few numerical results are presented to demonstrate the robustness and efficiency of the PSD solver.

Figures (5)

• Figure 1: The discrete $\ell^\infty$ and $\ell^2$ numerical errors vs. the iteration number, with a spatial resolution $N=256$. The numerical results are obtained by the proposed PSD iteration solver. The time step size and surface diffusion parameters are taken as: ${\Delta t}=0.01$, $\varepsilon=0.05$.
• Figure 2: Left: The number of iterations needed to obtain a machine precision for the PSD solver, in terms of $\varepsilon=0.01:0.01:0.1$, with a fixed ${\Delta t}=0.01$. Right: The number of iterations needed to obtain a machine precision for the PSD solver, in terms of ${\Delta t}=0.01:0.01:0.1$, with a fixed $\varepsilon=0.05$.
• Figure 3: The discrete $\ell^2$ and $\ell^\infty$ numerical errors versus temporal resolution $N_T$ for $N_T = 100:100:1000$, with a spatial resolution of $N=512$. The numerical results are obtained by the computation using the proposed PSD iteration solver to the numerical scheme \ref{['eqn:scheme']}. The surface diffusion parameter is taken to be $\varepsilon=0.5$, and the expansive parameter is set as $\theta_0 = 2$. The data lie roughly on curves $CN_T^{-1}$ for appropriate choices of $C$, confirming the full first order accuracy of the scheme.
• Figure 4: (Color online.) Snapshots of the phase variable at the indicated time instants over the domain $\Omega = [0,1]^2$, $\varepsilon = 0.005$, $\theta_0=3$, with a constant mobility ${\cal M} \equiv 1$.
• Figure 5: Log-log plot of the temporal evolution of the discrete energy for $\varepsilon=0.005$, $\theta_0=3$, with a constant mobility ${\cal M} \equiv 1$. The energy decreases similar to $a_e t^{b_e}$ until saturation. The red line represents the energy plot obtained by the simulations, while the straight blue line is obtained by least squares approximations to the energy data. The least squares fit is only taken for the linear part of the calculated data, and only up to $t=100$. The fitted line has the form $a_e t^{b_e}$, with $a_e = 0.01933$, $b_e=-0.3271$.

Theorems & Definitions (31)

• Proposition 2.1
• Proposition 2.2
• proof
• Theorem 2.3
• Theorem 2.4
• Lemma 2.5
• proof
• Lemma 2.6
• proof
• Lemma 3.1