Newton Iterations for Implicit Spatial Operators

An implicit method applied to a nonlinear PDE produces a nonlinear algebraic system per time step. The step cannot be computed by substitution; it requires an iteration whose convergence is not guaranteed by the stability of the time discretization alone. The iteration’s structure; what the Jacobian looks like, how it is assembled, and how its sparsity is exploited; determines whether an implicit scheme is practical on a real grid.

From the evolution equation to a residual

Write the semi-discrete system as $M\dot u = F(u)$ where $u \in \mathbb{R}^N$ holds the values at all grid points, $M$ is the mass matrix (the identity for a finite-difference discretization, a sparse mass matrix for finite elements), and $F$ includes the spatial operator. Backward Euler requires

$$ R(u^{n+1}) \equiv M\left(u^{n+1} - u^n\right) - \Delta t, F(u^{n+1}) = 0 $$

Solving $R(u) = 0$ for the unknown vector is the whole job of the step. Newton’s method linearizes around the current iterate $u^{(k)}$:

$$ J^{(k)},\delta u = -R(u^{(k)}), \qquad u^{(k+1)} = u^{(k)} + \delta u, \qquad J = M - \Delta t,J_F $$

where $J_F = \partial F/\partial u$ is the Jacobian of the spatial operator. For a problem in $d$ space dimensions with $n$ points per direction, $N = n^d$, and the dense cost of factorizing $J$ is $O(N^3)$, which for $n = 100$ in one dimension is $10^{12}$ operations; impossible. The sparsity is not an optimization; it is the only thing that makes the method feasible.

What the Jacobian’s sparsity pattern is

For a finite-difference discretization the pattern follows the stencil. A five-point Laplacian in two dimensions on a grid numbered row by row gives a matrix with five diagonals: the main diagonal, the two adjacent diagonals within a row, and two diagonals offset by $n_x$ for the row above and below. The number of nonzeros is $5N$, so the matrix has $O(N)$ entries against $N^2$ for a dense one, and a sparse factorization’s cost is governed by the fill-in that the elimination order creates rather than by $N$ alone.

The fill-in is the whole story of cost. Eliminating a variable couples its neighbors, and in two or three dimensions the neighbors form a graph whose elimination front grows. A banded ordering on an $n \times n$ grid gives a factorization cost of $O(n^4) = O(N^2)$ and storage $O(n^3)$, while a nested dissection ordering reduces the cost to about $O(N^{3/2})$ and the storage to about $O(N \log N)$. The specific scaling depends on the discretization and the ordering, but the qualitative point is stable: the ratio of implicit to explicit cost per step grows with dimension, which is what sets the point where a different method, such as an iterative solver, multigrid, or operator splitting, becomes the right one.

For a finite-volume discretization on an unstructured mesh, the pattern is not a fixed set of diagonals but a graph derived from cell adjacency. The assembly is by flux, one contribution per face, and the matrix is assembled as a sparse structure rather than filled by index arithmetic.

Assembling the Jacobian without inventing entries

Two methods produce $J_F$. The analytic method differentiates the flux function and the stencil by hand. For a flux $G(u, \nabla u)$ the derivative has a term $\partial G/\partial u$ acting on point values and a term $\partial G/\partial(\nabla u)$ acting on differences, so each face contributes a small block to the matrix. This is exact, exploits the flux’s structure, and reproduces the sparsity pattern in the assembly code.

The numerical method perturbs each unknown:

/* Column-wise finite-difference Jacobian, dense storage. */
double eps = 1e-7;                               /* relative perturbation */
for (int j = 0; j < N; j++) {
    double saved = u[j];
    double h = eps * (fabs(saved) > 1.0 ? fabs(saved) : 1.0);
    u[j] = saved + h;
    compute_F(u, Fp);
    u[j] = saved;
    for (int i = 0; i < N; i++)
        J[i][j] = (Fp[i] - F[i]) / h;            /* one column at a time */
}

This costs $N + 1$ evaluations of $F$ for a dense matrix, and for a large problem that is far more expensive than the factorization it feeds. The sparsity helps: since column $j$ can only have nonzeros in its stencil neighbors, a coloring of $f$‘s dependency graph permits several columns to be perturbed simultaneously, and each evaluation then fills a whole group of columns. For nearest-neighbor stencils the number of colors is small and the assembly cost becomes a modest multiple of $n$ evaluations rather than $N$. The correctness requirement is that no two columns of the same color share a row, otherwise their contributions are summed in a single evaluation and cannot be separated.

The perturbation size is a genuine trade-off. Too large and the difference quotient is not the derivative; too small and the subtraction loses digits. A relative perturbation around the square root of machine epsilon is the conventional compromise, and the resulting entries carry a few digits of relative error, which is harmless for a Newton iteration whose convergence is quadratic only up to the accuracy of the linearization.

Convergence and what breaks it

Newton’s method converges quadratically near a solution with a nonsingular Jacobian: the error satisfies $e_{k+1} \le c,e_k^2$. Far from it, the guarantee is local at best. Three failures dominate in practice.

The first is a bad initial guess. The previous time step’s solution is the natural starting value, and it is accurate when the step is small; when the step is large, as it must be in a stiff problem, the guess is far from the new solution and the iteration may diverge. A damped Newton, which shortens the update when the residual norm does not decrease, restores global convergence for many problems but not all.

The second is a singular or ill-conditioned Jacobian. As $\Delta t \to \infty$ in a steady state problem, $J \approx -\Delta t J_F$, and if $J_F$ has a zero eigenvalue, as a conserved quantity with no dissipation does, then the system is singular. The time term $M$ regularizes this for finite $\Delta t$, but the conditioning degrades as $\Delta t$ grows. The mechanism of the singularity is the same one that makes the steady problem underdetermined.

The third is a nonsmooth right-hand side. A limiter, a clipped reaction rate, or a phase change makes $F$ continuous but not differentiable. The accumulated Jacobian is the derivative of one branch, and the iteration can cycle. The standard remedies are to smooth the nonsmooth term, to use a semismooth variant, or to accept indifferentiability and rely on damping.

Inexact Newton and the two tolerances

The Jacobian solve does not need to be exact. The inexact Newton framework replaces the exact update with one that satisfies

$$ |J\delta u + R| \le \eta_k |R| $$

for a forcing term $\eta_k$. With $\eta_k$ chosen proportional to the residual’s size, the outer convergence stays superlinear, and the total work is often much less than solving each system to precision. This matters because the inner solve, by an iterative linear method, is itself the dominant cost.

The two tolerances on the step are then separate and should not be conflated. The Newton tolerance decides when the nonlinear system is solved; the step’s accuracy relative to the time discretization’s error is what actually matters. Solving Newton to machine precision when the time discretization error is $10^{-5}$ is wasted work. A common and defensible choice is a Newton tolerance a small fraction of the time discretization’s local error estimate, with a floor at the roundoff level. Where the Newton solve’s accuracy interacts with the step-size controller is worth stating explicitly: an inaccurate Newton solve leaves a residual in the discrete equation, the error estimate computed from the step then reflects that residual as well as the truncation error, and a controller comparing that estimate to its target responds by cutting the step. The cut is the right response only if the inaccuracy is genuine, which is why the Newton tolerance should sit below the tolerance handed to the controller by a comfortable margin. The same separation appears in the spatial operator of an implicit diffusion scheme, where the linear solve is the whole step and the only tolerance is the factorization’s; see explicit and implicit treatment of diffusion for the linear case and the structure of the matrix it factors.

The forcing term needs a floor for the same reason the step size does. As $|R|$ shrinks toward the roundoff level of the residual evaluation, the condition $|J\delta u + R| \le \eta_k |R|$ becomes unachievable, and an iterative inner solver will either fail to converge or spend its whole iteration budget on noise. A forcing term that is a fixed fraction rather than a power of the residual, or an absolute floor on the requested reduction, avoids the wasted work.

There is a limit on how much the inner solve can be relaxed before the outer iteration loses its superlinear rate. Choosing $\eta_k$ proportional to the previous residual’s size is the standard compromise: the forcing term shrinks as the outer iterate improves, so the last iterations are solved accurately and the early ones are not. The choice of proportionality constant decides where on the spectrum between cheap-but-slow and expensive-but-fast the method sits, and it is the parameter that a stalled nonlinear solve most often turns out to have set badly.

Reusing the Jacobian across steps

A Newton Jacobian at one time step is close to the Jacobian at the next when the step is small and the solution moves smoothly. Reusing the factorization saves the most expensive part of the step, and a simple-minded implementation can pay for the reuse by failing to converge. The usual compromise is a reuse limit: keep the factorization for a handful of steps, then rebuild. The modified Newton method converges linearly with iteration matrix whose spectral radius is the perturbation in the Jacobian, so the reuse limit is chosen so that the linear rate does not exceed the number of Newton iterations budgeted. A direction-aware version, in which the convergence rate is used to estimate the effective spectral radius, allows reuse to extend when the problem is behaving.

For a linear PDE the Jacobian is constant, and the entire discussion collapses: one factorization serves every step. This is why linear diffusion and advection solvers are much cheaper per step than nonlinear ones, and why the linear case in discretization and stability needs only the solve.