Transform before linearizing: robust Newton methods for singular $p$-Laplace and $p$-Stokes equations

1 Introduction

• Newton’s fast convergence is confined to iterates near the solution, while Picard is more globally robust but converges only linearly, with its contraction worsening as p approaches 1.For Newton, iteration counts typically grow under mesh refinement and depend severely on p.Source: §1 · p. 2

1.1 Approach and related literature

• The method lifts the equation with an auxiliary variable, transforms the lifted system before linearizing, then eliminates that variable by static condensation.For p-Laplace, the lifting variable is the flux; elimination keeps the cost per iteration comparable to standard Newton or Picard.Source: §1.1 · pp. 2–3
• The authors position their method among lifted and primal–dual approaches, while distinguishing it from preconditioning that targets locally difficult subdomains.They also connect stress lifting in related PDE solvers to reformulating constitutive relations with auxiliary variables.Source: §1.1 · p. 3

1.2 Main contributions and limitations

• The method yields iteration counts largely insensitive to p and mesh refinement, with up to an order-of-magnitude fewer linear solves than standard Newton for p-Stokes and time-dependent p-Navier–Stokes near p = 1.The reported p-Laplace experiments also cover two and three dimensions.Source: §1.2 · p. 4
• The method is most advantageous for strongly nonlinear, potentially singular constitutive laws, and its convergence analysis covers only discretized, regularized p-Laplace equations.The analysis constants deteriorate as p approaches 1, and the experiments do not study preconditioned iterative Krylov solvers.Source: §1.2 · p. 4

2 The p-Laplace equation

• For homogeneous boundary data, the weak p-Laplace solution is the unique minimizer of a strictly convex energy in W^1,p.Source: §2 · p. 4

2.1 Dual problem and extremality conditions

• The flux lifting makes the PDE linear in the primal variable and confines nonlinearity to a pointwise constitutive relation, with equivalent forms yielding different Newton linearizations.Source: §2 · pp. 4–5, §2.1 · p. 5

2.2 Regularization and Newton methods

• Lifting alone recovers standard Newton, whereas nonlinear rescaling before linearization creates a distinct step using the lagged dual direction instead of the primal direction.Source: §2.2 · pp. 6–7

2.3 Modifications to ensure well-posedness

• The pointwise bound |wλ| ≤ β/(2-p), with 2-p < β < 1, keeps diffusion uniformly positive definite and becomes inactive near the solution.Source: §2.3 · p. 7, §2.4 · p. 8

2.4 Summary of lifted Newton algorithm

• Projection acts only in constructing the diffusion tensor, and λ0 = 0 makes the first step Picard; each iteration requires only one linear solve.Source: §2.4 · pp. 8–9

3.1 Global convergence for discretized method

• The discretized lifted iteration converges globally from arbitrary initializations, with residuals tending to zero and iterates reaching the unique minimizer of the regularized energy.The proof uses a finite-dimensional energy sublevel set and line search; its constants may depend on the discretization.Source: §3.1 · p. 9
• The global-convergence coercivity bound is not mesh-independent: it relies on a finite-dimensional gradient bound unavailable from continuous energy sublevels alone.The analysis constants also deteriorate as ε→0 and p→1, at the same rate as for standard Newton.Source: §3.1 · p. 11, §3.2 · p. 13

3.2 Local convergence for discretized problem

• Once the combined primal error and lifting inconsistency are small, the method takes unit steps and converges quadratically.The lifting update cancels the first-order change in the constitutive relation up to Jacobian error, reducing inconsistency quadratically.Source: §3.2 · pp. 11–12

3.3 Robustness away from the solution

• In regions with |∇u| ∼ ε, the standard Newton direction can vary sharply between iterations, while a lagged lifting direction stays close to Picard and filters those fluctuations.This effect is confined to iterations before the Newton convergence basin is reached.Source: §3.3 · p. 13
• For the one-dimensional p-Laplace example at p=1.01, standard Newton zig-zags and takes 27 iterations, while lifted Newton converges in 4.With a small lagged λ, the lifted update sets plateau gradients to zero up to O(p−1) in one step.Source: §3.3 · pp. 13–14, Figure 1 · p. 14

4.1 Implementation and comparisons

• The lifting variable is updated algebraically at quadrature points, so it is not represented as a finite-element function and requires no inf-sup condition.The compared methods include Picard, standard Newton, lifted Newton, and Picard–Newton; iteration counts correspond to linear solves, with each Picard–Newton step counted twice.Source: §4.1 · pp. 14–15

4.2 Two-dimensional p-Laplace equation example

• As p approaches 1, two-dimensional iteration counts rise steeply for all methods except lifted Newton, whose counts grow only mildly across all three grids.Picard fails to reach tolerance for p = 1.01, and lifted Newton is less sensitive to mesh refinement; standard Newton accepts many strongly damped steps while lifted Newton soon accepts full steps.Source: §4.2 · p. 16, Figure 2 · p. 16, Table 2 · p. 17
• For p near 1, decreasing ε generally increases iteration counts, and at ε = 10^-10 with p = 1.01 both methods fail to reach the prescribed tolerance.The authors attribute these failures likely to roundoff errors at the extreme parameter values and identify ε-regularization as necessary for solver robustness.Source: §4.2 · p. 17, Table 3 · p. 17

4.3 Three-dimensional p-Laplace equation example

• For p = 1.01, lifted Newton needs 21 linear solves with CG1 and 24 with CG2, versus 257 and 242 for standard Newton on the three-dimensional mesh.Picard does not meet tolerance within 500 iterations for either discretization; Picard–Newton takes 50 and 44 linear solves, respectively.Source: §4.3 · p. 18, Table 4 · p. 19

• For 1 < p < 2, p-Stokes describes shear-thinning power-law flow, and its near-singular viscosity can make Newton slow and fragile as p approaches 1.Unlike the p-Laplace problem, its unknown is a divergence-free vector field coupled to pressure.Source: §5 · p. 19
• The method lifts the regularized stress relation with the deviatoric stress tensor, then linearizes and eliminates its increment to obtain a velocity–pressure system.The lifting variable influences the velocity increment through the diffusion tensor M.Source: §5 · pp. 19–20
• Symmetrizing the diffusion tensor preserves its quadratic form, while a pointwise feasibility bound on the stress direction makes the velocity block uniformly elliptic.The resulting system is a symmetric saddle-point system, allowing standard Stokes solvers and preconditioners.Source: §5 · p. 20

6 Numerical results: p-Stokes

• The p-Stokes experiments use Taylor–Hood elements and compare standard and lifted Newton under mesh refinement and in time-dependent flow.The time-dependent case uses implicit Euler, Nitsche boundary conditions, and continuous interior penalty stabilization.Source: §6 · p. 20, §6.2 · p. 21

6.1 Lid-driven cavity flow

• For p = 1.01 on the finest cavity mesh, lifted Newton takes 22 iterations versus 285 for standard Newton.At p = 1.5, lifted Newton takes four or five iterations on all three meshes, while standard Newton rises from 19 to 47; both methods show mild mesh dependence at p = 1.01.Source: §6.1 · p. 21

6.2 Flow around a cylinder for the p-Navier–Stokes equations

• At p = 1.01 in the cylinder-flow simulation, lifted Newton reduces the mean iterations per time step from 9.86 to 4.12 and the maximum from 140 to 24.Across 3,392 time steps, the reported linear solves fall from roughly 33,400 for standard Newton to 14,000 for lifted Newton; for t > 4, both methods require similar iterations.Source: §6.2 · p. 22

Continuation

• No reported experiments use continuation; the authors say it is unnecessary because lifted iteration counts are largely insensitive to p and mesh refinement.Continuation would require solving a sequence of auxiliary problems, with total cost across intermediate solves.Source: §7 · pp. 22–23

Krylov solvers

• For larger three-dimensional problems, the lifted tangent’s spectral equivalence to the Picard operator makes Picard-equivalent preconditioners candidates for iterative solves.The reported linear systems were solved directly.Source: §7 · p. 23

Symmetrized versus non-symmetrized lifted Newton operator

• The symmetrized and nonsymmetrized lifted operators have identical quadratic forms, coercivity bounds, descent properties, and global convergence guarantees.For p = 1.01 and N = 25 in the two-dimensional p-Laplace example, symmetrization required 23 iterations versus 40 without it.Source: §7 · p. 23

The case p > 2

• For p > 2, the proposed transformation is unsuitable; the authors suggest lifting the dual constitutive relation instead and leave systematic study for future work.They cite degeneracy of the diffusion coefficient and ill-conditioning as distinct difficulties in this regime.Source: §7 · pp. 23–24

Open questions

• The convergence bounds deteriorate as p approaches 1, and local analysis quantifying the observed robustness, infinite-dimensional convergence, and the limit ε → 0 remain open.The authors note that the infinite-dimensional challenge is shared with standard Newton.Source: §7 · p. 24

A Alternative projection of the lifting variable

• The alternative projects the lifting iterate itself, making the projection in the diffusion-tensor definition unnecessary once the feasibility bound holds for that iterate.The update uses the projected value pointwise.Source: A Alternative projection of the… · p. 26
• The two variants differ in the initial phase, but the analysis covers both and the flux-convergence conclusion is unchanged.The pointwise scaling projection is nonexpansive and preserves the stated recursion.Source: A Alternative projection of the… · p. 26
• The two variants yield almost identical iteration counts in practice, as summarized for the two-dimensional p-Laplace example on 10,000 triangles.The table caption identifies the comparison as lifted Newton versus the alternative projection variant.Source: A Alternative projection of the… · p. 26, Table 5 · p. 27

B One-dimensional finite difference example

• Standard Newton uses a nearly flat tangent slope when p is close to 1 and overshoots by a factor of 1/(p − 1), whereas lifted Newton uses a secant slope and gets close to the solution quickly.The comparison concerns the conversion of the shared flux increment into a gradient increment.Source: B One-dimensional finite difference… · p. 27
• The finite-difference problem uses homogeneous Dirichlet conditions on (0, 1), with uniformly spaced nodes and discrete gradient y = Bu.The forcing is chosen so the discrete solution is constant at the interior degrees of freedom, producing a plateau except on the two boundary intervals.Source: B One-dimensional finite difference… · p. 27
• Both methods solve a linear system with the same residual and differ in its diagonal coefficients; the lifted version uses a feasibility projection for the lifting variable.The corresponding lifted diffusion tensor is described as the one-dimensional counterpart of the tensor in equation (27).Source: B One-dimensional finite difference… · p. 27

B.1 Expansions for small values of p −1

• With τ = p − 1 ≪ 1 and gradients bounded away from zero and infinity, the constitutive function is nearly the sign function: its tangent slope is order τ, while its secant slope is order one.The analysis assumes |y_i| ≫ ε, so |y_i|_ε ≈ |y_i|.Source: §B.1 · p. 28
• The lifted method uses the secant slope when π_u(λ_i) = O(τ), while setting λ = g(y) recovers the tangent slope.This explains why the choice of lifting variable changes the gradient update.Source: §B.1 · p. 28
• In one dimension, the flux increment is determined by the residual up to one constant fixed by the homogeneous boundary conditions, yielding an explicit solution formula for any positive diagonal D.The flux increment is defined as ẑ = D B û, and the formula follows from the residual relation and boundary conditions.Source: §B.1 · p. 28

B.2 The flux increment is the same for both methods

• In the flat interior, hρ_i ≈ −g(y_i) ≈ −sign(y_i) for i = 2, …, n, while the boundary residuals satisfy hρ_1 = 0 and hρ_{n+1} = O(τ).These estimates hold close to the solution.Source: §B.2 · p. 28
• The estimate uses δ := max_i |u_i − 1| as the iterate’s distance from the solution and assumes δ ≪ 1.The derivation also uses the stated interior and boundary estimates for d^{-1}_i.Source: §B.2 · pp. 28–29
• For i = 2, …, n, both methods produce the flux increment ẑ_i ≈ −g(y_i) ≈ −sign(y_i).The resulting linearized flux vanishes to leading order in the flat region.Source: §B.2 · p. 29

B.3 Comparison of the primal increments

• A unit standard Newton step changes each flat-region gradient’s sign and increases its magnitude by approximately 1/τ, producing zig-zagging.The residual then changes sign and the next step reverses the motion; a line search must compromise between the flat region and steep boundary intervals.Source: §B.3 · p. 29
• A unit lifted Newton step sets flat-region gradients to zero up to O(τ).The conversion from flux increment to gradient increment uses the secant slope rather than the nearly flat function’s tangent slope.Source: §B.3 · p. 29

B.4 The lifting variable in the flat region

• The lifting-variable update is exactly the linearized flux after the step, so its flat-region components are O(τ) for every prior lifting variable.This follows from the flat-region flux increment estimate, which does not depend on D.Source: §B.4 · p. 29
• With λ_0 = 0, the first lifted Newton step coincides with a Picard step and equals τ times the first standard Newton step.The flat-region estimate repeats after each step.Source: §B.4 · p. 29
• At the boundary intervals, λ_1 ≈ 1 and λ_{n+1} ≈ −1 after the first step, so lifted Newton reduces to standard Newton, which converges quadratically in the steep region.Source: §B.4 · p. 29
• The lagged dual direction filters the extreme sensitivity of w_ε where the gradient is small.In the one-dimensional model, the secant slope makes the lifted step land on the flat solution up to O(τ), unlike standard Newton’s tangent-slope conversion.Source: §B.4 · p. 29