Projection-free approximation of flows of harmonic maps with quadratic constraint accuracy and variable step sizes

1. Problem Statement

The paper addresses the numerical approximation of harmonic map flows into spheres. These flows are gradient flows of the Dirichlet energy functional:
$$I[u] = \frac{1}{2} \int_{\Omega} |\nabla u|^2 \, dx$$
subject to the unit-length constraint $|u(x)| = 1$ almost everywhere in the domain $\Omega$. The evolution is governed by the harmonic map heat flow equation:
$$\partial_t u - \Delta u = |\nabla u|^2 u$$
with initial and boundary conditions.

Key Challenges:

• Constraint Preservation: Standard numerical methods often fail to preserve the unit-length constraint $|u|=1$ exactly at the discrete level, leading to constraint violations that can accumulate.
• Projection Costs: Methods that enforce the constraint exactly (via projection) often require solving nonlinear systems or impose strict restrictions on mesh types and step sizes.
• Accuracy vs. Stability: Existing projection-free linearly implicit methods (like the standard Euler scheme) are unconditionally energy stable but typically only achieve first-order accuracy regarding the constraint violation.
• Variable Step Sizes: Many high-order methods (like BDF2) struggle with unconditional stability when variable step sizes are introduced, limiting their ability to adapt to singularities or accelerate convergence.

2. Methodology

The authors propose a projection-free linearly implicit $(\theta, \mu)$-method based on Alouges' weak formulation, where the test space is restricted to be orthogonal to the solution ($v \cdot u = 0$).

The $(\theta, \mu)$-Method

The method iteratively computes approximations $u^n$ using a time step $\tau_n$. For $n \geq 2$, it seeks a discrete time derivative $d_t u^n$ satisfying a linearized unit-length condition $d_t u^n \cdot \tilde{u}^{n-1+\mu} = 0$, where $\tilde{u}$ is an extrapolated value. The variational formulation is:
$$(d_t u^n, v)_* + (\nabla [u^{n-1} + \theta \tau_n d_t u^n], \nabla v) = 0$$
for all test functions $v$ orthogonal to the extrapolated value $\tilde{u}^{n-1+\mu} = u^{n-1} + \mu \tau_n d_t u^{n-1}$.

• Parameters:
  • $\theta \in (0, 1]$: Controls energy dissipation.
  • $\mu \in [0, 1]$: Controls the point of orthogonality for the constraint.
• Special Cases:
  • $\theta=1, \mu=0$: Reduces to the standard linearly implicit Euler method (first-order constraint accuracy).
  • $\theta=1/2, \mu=1/2$: The Linearly Implicit Midpoint Method (proposed as the optimal choice).
  • $\theta=1, \mu=1/2$: A Modified Implicit Euler method (uses Euler's variational form but midpoint's constraint orthogonality).

Variable Step Size Extension

The authors extend the method to variable step sizes ($\tau_n$). They derive a discrete chain rule and energy identities that hold even when $\tau_n$ varies. This allows for:

1. Acceleration: Increasing step sizes as the solution approaches a stationary state.
2. Adaptivity: Decreasing step sizes near singularities (e.g., blow-up or non-smooth initial data) to maintain accuracy.

3. Key Contributions

1. Quadratic Constraint Accuracy:
   The paper proves that the Midpoint Method ($\theta=1/2, \mu=1/2$) achieves second-order accuracy in the constraint violation ($|u|^2 - 1 = O(\tau^2)$) under a sharp discrete regularity condition. This is a significant improvement over the standard Euler method, which is only $O(\tau)$.

2. Unconditional Energy Stability:
   The proposed method is proven to be unconditionally energy stable (energy decreasing) for any step size, regardless of whether the step size is constant or variable. This is a crucial property for long-time integration.

3. Variable Step Size Stability:
   Unlike many high-order methods (like BDF2) where variable step sizes can destroy stability, the $(\theta, \mu)$-method retains its unconditional stability and constraint properties with variable steps. The authors provide a specific update rule ($\tau_{n+1} = \tau_n \sqrt{1 + c\tau_n}$) to balance error terms and accelerate convergence.

4. Comparison with BDF2:
   The paper provides a rigorous theoretical and numerical comparison with the projection-free two-step BDF method. While BDF2 also offers second-order accuracy, it requires a complex "G-stability" analysis and is more sensitive to variable step sizes. The $(\theta, \mu)$-method offers a simpler implementation with comparable or superior performance in specific metrics.

4. Numerical Results

The authors conducted extensive numerical experiments in 2D using finite element discretization:

• Constant Step Size (H1 and L2 Gradient Flows):

  • The Midpoint Method and Modified Euler consistently demonstrated second-order convergence in the $L^1$-norm of the constraint violation error, whereas the standard Euler showed only first-order convergence.
  • The Midpoint method was slightly more accurate than the BDF2 method (constraint error was $\approx 1/3$ of BDF2's).
  • Energy approximation errors stabilized at a level determined by spatial discretization, confirming the temporal accuracy was not the limiting factor.

• Singularities:

  • In tests with initial data leading to finite-time singularities, the method maintained stability. The convergence rate dropped to superlinear (but not quadratic) due to the lack of solution regularity, consistent with theoretical expectations.

• Variable Step Sizes:

  • Acceleration: Using a prescribed increasing step size rule significantly reduced the constraint violation error for the same computational cost (number of iterations) compared to constant step sizes in the H1-gradient flow case.
  • Adaptivity: An adaptive step size control mechanism (reducing $\tau$ when the update norm increases) successfully handled singularities. It reduced the total number of iterations by 20% compared to constant step sizes while maintaining the same accuracy, by automatically refining the mesh in time during rapid dynamics and coarsening it during slow evolution.

5. Significance

This work is significant for the numerical analysis of geometric PDEs because:

• It bridges the gap between efficiency (linearly implicit, no nonlinear solves) and accuracy (second-order constraint preservation).
• It solves a major limitation of projection-free methods: the inability to use variable step sizes without sacrificing stability. This is critical for practical applications where solutions may evolve rapidly (singularities) or slowly (convergence to equilibrium).
• It provides a robust alternative to BDF methods, offering simpler implementation (no need for G-stability matrices) and superior performance in constraint preservation for harmonic map flows.
• The methodology is generalizable to other geometrically constrained problems, such as the Landau-Lifshitz-Gilbert (LLG) equation in micromagnetics and liquid crystal models.