Nonlinear Solver Methods: Newton, Broyden, and Quasi-Newton for PDEs

When you discretize a nonlinear partial differential equation, you get a nonlinear algebraic system. Solving it requires choosing the right solver method. Newton’s method converges quadratically but demands a full Jacobian. Broyden’s method builds the Jacobian cheaper but can destabilize. The Jacobian-free Newton-Krylov (JFNK) approach avoids the Jacobian entirely, relying only on Jacobian-vector products with physics-based preconditioning to keep convergence tractable.

The choice is not theoretical. It determines whether your simulation converges at all, how much memory it consumes, and whether it finishes in minutes or days.

Link to existing article: Time Integration Methods for PDE Solvers covers the broader landscape of explicit versus implicit schemes and operator splitting. Link to existing article: Implicit vs Explicit Methods discusses stability regions, accuracy criteria, and when implicit schemes become necessary. This article focuses specifically on the nonlinear solvers that sit inside every implicit scheme.

When Implicit Methods Demand a Nonlinear Solver

Implicit time integration requires solving a nonlinear system at every time step. Consider the backward Euler discretization of a generic PDE:

F(uⁿ⁺¹) = 0

The unknown uⁿ⁺¹ appears inside the operator F through diffusion terms, reaction kinetics, or coupled multi-physics equations. After spatial discretization (finite volume, finite element, or finite difference), the system becomes:

F(x) = 0, where x ∈ ℝⁿ

The dimension n is the number of degrees of freedom — typically the product of grid cells, physical fields, and time-step iterations. In a 3D phase-field simulation with 500³ cells and two coupled fields, n can easily exceed 10⁸.

The challenge is that F is nonlinear. You cannot write x = F⁻¹(0). Instead, you need an iterative algorithm that constructs a sequence converging to the root.

Newton’s Method: The Baseline

Newton’s method is the reference against which all other nonlinear solvers are measured. For the system F(x) = 0, it constructs the Jacobian matrix:

J = ∂F/∂x

And iterates:

xₖ₊₁ = xₖ − J⁻¹·F(xₖ)

The key property is quadratic convergence. If the Jacobian is Lipschitz continuous and the initial guess is within the convergence basin, the error satisfies:

‖xₖ₊₁ − x*‖ ≤ C · ‖xₖ − x*‖²

This means the number of correct digits roughly doubles with each iteration. For well-behaved PDE systems, Newton’s method often converges in 3–5 iterations.

Full Jacobian: Analytical vs Finite Difference

Constructing the Jacobian is the most expensive part of Newton’s method. The options are:

Analytical Jacobian. You derive ∂F/∂x by hand or through automatic differentiation. In libraries like FEniCS, the variational formulation automatically generates the Jacobian through symbolic differentiation. The result is exact (up to floating-point precision).

Finite-difference Jacobian. You perturb each column of J by ε and re-evaluate F:

Jᵢ ≈ (F(x + ε·eᵢ) − F(x)) / ε

This requires n additional function evaluations per Newton iteration. For large n, this dominates the cost.

Block Jacobian. In coupled systems with k fields, the Jacobian is block-structured:

J = [ ∂F₁/∂x₁  ∂F₁/∂x₂  ... ]
    [ ∂F₂/∂x₁  ∂F₂/∂x₂  ... ]
    [          ...       ... ]

If you only need diagonal blocks ∂Fᵢ/∂xᵢ, you can evaluate them independently — a technique called diagonal freezing or fully coupled vs segregated solvers.

Global Convergence: Line Search and Trust Regions

Pure Newton’s method is a local method. If the initial guess is far from the solution, the iterations can diverge even though quadratic convergence holds near the root.

Line search modifies the Newton step by scaling it:

xₖ₊₁ = xₖ − α·J⁻¹·F(xₖ)

where α ∈ (0, 1] is chosen to reduce a merit function. The classic approach minimizes ‖F(xₖ − α·J⁻¹·F(xₖ))‖.

Trust regions define a radius Δ within which the quadratic model is trusted:

min q(s) = F(xₖ) + J·s + ½·sᵀ·H·s
subject to ‖s‖ ≤ Δ

When the reduction in the merit function is insufficient, the radius shrinks. When it is sufficient, the radius grows.

PETSc’s SNES (SNU Nonlinear Equation Solver) package provides both line-search and trust-region globalization. The choice affects the total iteration count more than the per-iteration cost.

Practical Newton: When It Works and When It Doesn’t

Newton’s method excels when:

• The Jacobian is available in exact form (symbolic, AD, or block-analytical)
• The nonlinearities are mild (diffusion, mild reaction terms)
• Memory is sufficient for the full Jacobian (or its sparse factorization)

Newton’s method struggles when:

• n exceeds 10⁶ and assembling the full Jacobian is prohibitive
• The function F is opaque (black-box simulation code)
• Each function evaluation already dominates the runtime

For large-scale PDEs, the full Jacobian approach is often too expensive. This is where quasi-Newton methods enter.

Broyden’s Quasi-Newton Methods

Broyden’s method belongs to the broader class of quasi-Newton methods. Instead of computing J = ∂F/∂x, it builds an approximation Bₖ ≈ J through rank-1 updates to the secant condition:

Bₖ₊₁ · (xₖ₊₁ − xₖ) = F(xₖ₊₁) − F(xₖ)

The secant condition is one equation for n unknowns in Bₖ₊₁. Broyden’s formulation chooses the update that minimizes the Frobenius norm change:

min ‖Bₖ₊₁ − Bₖ‖_F
subject to Bₖ₊₁·(xₖ₊₁ − xₖ) = F(xₖ₊₁) − F(xₖ)

Good vs Bad Broyden

Broyden’s original method (sometimes called “good Broyden”) updates Bₖ using the secant equation. The inverse variant (“bad Broyden”) updates Bₖ⁻¹ directly:

Bₖ₊₁⁻¹ = Bₖ⁻¹ + (δy − Bₖ⁻¹·δx)·δxᵀ / (δxᵀ·δx)

where δx = xₖ₊₁ − xₖ and δy = F(xₖ₊₁) − F(xₖ).

The inverse update is what Python’s scipy.optimize.broyden1 uses. It avoids inverting the Jacobian approximation at each step — a significant saving for large problems.

The Condition-Number Instability

The most important practical limitation of Broyden’s method is that its Jacobian approximation can become badly conditioned. The rank-1 update modifies only one direction, while the rest of the matrix drifts. Over many iterations, Bₖ accumulates outdated curvature information.

The condition number κ(Bₖ) = ‖Bₖ‖ · ‖Bₖ⁻¹‖ can grow without bound. A condition number exceeding 10⁸ causes Krylov solvers like GMRES to stall, even if the underlying problem is well-conditioned.

This is the concrete reason to prefer Newton or JFNK over Broyden in production simulations. For mildly nonlinear problems with few degrees of freedom per field (e.g., 1D diffusion, 2D steady-state), Broyden’s method can be competitive. For 3D systems with 10⁶+ degrees of freedom, the condition drift makes it unreliable.

Limited-Memory Broyden (L-BFGS)

For large-scale optimization, L-BFGS stores the last m step pairs (δxᵢ, δyᵢ) and implicitly constructs the inverse Jacobian through a recursion. It avoids storing the full n × n matrix.

However, L-BFGS is designed for optimization (minimizing a scalar objective), not root-finding (solving F(x) = 0). For PDE root-finding, the limited-memory Broyden update used in solvers like PETSc’s BFGS solver (not L-BFGS) is more common.

Jacobian-Free Newton-Krylov: The Modern Large-Scale Standard

The Jacobian-free Newton-Krylov (JFNK) method addresses the memory bottleneck of full Jacobian Newton by never assembling the Jacobian. Instead, it uses a Krylov subspace method (typically GMRES) to solve the Newton system:

J · Δx = −F(x)

GMRES only requires matrix-vector products J·v. These are computed through the Jacobian-vector product:

J · v ≈ (F(x + ε·v) − F(x)) / ε

No explicit Jacobian is assembled. Only F and the perturbation vector v are needed.

Why Preconditioning Is Mandatory

JFNK with a bare Jacobian-vector product is rarely practical. The Newton system can have a large condition number, especially when the Jacobian has widely varying eigenvalues from different physical processes.

Physics-based preconditioning (PBP) addresses this by factoring out known operators. Consider a convection-diffusion-reaction system:

J = J_convection + J_diffusion + J_reaction

PBP replaces the preconditioner with:

M = J_diffusion

The diffusion operator usually dominates the spectrum and is often block-tridiagonal or tridiagonal per field. Its inverse can be computed efficiently through sparse direct factorization or ILU, while the convection and reaction terms are treated explicitly.

The preconditioned system becomes:

M⁻¹ · J · Δx = −M⁻¹ · F(x)

The condition number of M⁻¹·J is dramatically reduced because the diffusion operator — the source of the worst conditioning — has been factored into M.

Without preconditioning, JFNK can require hundreds or thousands of GMRES iterations. With PBP, it often converges in 10–50 iterations, matching or exceeding full-Jacobian Newton.

The Newton-Krylov Pattern in Practice

1. Construct the Newton residual r = F(xₖ)
2. Solve J · Δx = −r using GMRES
3. Jacobian-vector products are computed through finite-difference perturbations of F
4. Preconditioner M is applied at each GMRES iteration
5. Update xₖ₊₁ = xₖ + Δx
6. Check convergence on ‖r‖

This pattern appears in every major PDE framework: FEniCS/DOLFIN, PETSc/SNES, FiPy, and NonlinearSolve.jl. The difference lies in how the preconditioner is constructed and how the Jacobian-vector product is computed.

The Polyalgorithm Pattern: How Modern Solvers Actually Deploy

Production solvers rarely rely on a single method. Instead, they use a polyalgorithm that adapts to problem difficulty:

1. Start with fast Broyden (cheap per iteration, no Jacobian)
2. If convergence stalls, fall back to full Newton (robust quadratic convergence)
3. If Newton also fails, fall back to TrustRegion (globally convergent)

This is the pattern used by PETSc’s SNES, SciML’s DifferentialEquations.jl, and MOOSE’s nonlinear solver.

Eisenstat-Walker Forcing

The Eisenstat-Walker strategy (sometimes called the “forcing strategy”) adapts the inner tolerance of the Krylov solver relative to the Newton residual:

‖J · Δx + r‖ ≤ ηₖ · ‖r‖

where ηₖ relaxes as r shrinks. Early Newton iterations use a loose inner tolerance (fewer Krylov steps). As the Newton residual decreases, the inner tolerance tightens (more Krylov accuracy).

This prevents wasted work solving the linear system to high accuracy when the Newton iterate is still far from the root.

Decision Framework: When to Choose Which Method

The following table summarizes the practical trade-offs:

Method Comparison

Choosing by Problem Size

Small problems (n < 10⁴): Full Newton with analytical or symbolic Jacobian. The cost is dominated by matrix assembly, so the quadratic convergence outweighs the setup cost.

Medium problems (10⁴ < n < 10⁶): Broyden can be competitive if the nonlinearity is mild. Newton with sparse direct factorization is more robust but costs more per iteration.

Large problems (n > 10⁶): JFNK with physics-based preconditioning is the standard choice. Full Jacobian assembly is prohibitive, and Broyden’s condition drift is unreliable.

Choosing by Nonlinearity Severity

Mild nonlinearity (diffusion-dominated, weak reaction terms): Newton or Broyden both work well. The quadratic convergence of Newton makes it attractive for small systems.

Strong nonlinearity (coupled phase-field, reaction-diffusion with stiff kinetics): Newton with globalization (line search or trust region) is essential. Broyden’s approximation can mislead the search direction.

Opaque nonlinearity (black-box simulation, unknown F): JFNK or Broyden, since they do not require explicit Jacobian construction.

Choosing by Memory Constraints

Limited memory (single-node, < 64 GB): JFNK requires only O(n) storage for vectors. Newton requires O(nnz(J)) for sparse Jacobian storage plus O(nnz(L)) for factorization.

Ample memory (cluster node, > 256 GB): Full Newton becomes feasible for problems with n up to 10⁷, provided the Jacobian is sparse.

Code Examples: How These Methods Look in Practice

SciPy: Broyden’s Method

SciPy provides broyden1 for root-finding with Broyden’s quasi-Newton update:

from scipy.optimize import broyden1, root

def f(x):
    """System F(x) = 0"""
    return [x[0]**2 + x[1] - 1,
            x[0] + x[1]**2 - 2]

result = broyden1(f, [0.5, 0.5])
print(f"Solution: {result.x}")

The broyden1 function builds the inverse Jacobian approximation through rank-1 updates. For small systems, it avoids full Jacobian assembly. For large systems, the condition drift can cause convergence failure.

The root interface also supports 'hybr' (MINPACK’s Levenberg-Trust-Region) and 'lm' (Levenberg-Marquardt), both of which can handle structured least-squares problems better than Broyden.

FEniCS: Newton with Symbolic Jacobian

FEniCS generates the Jacobian symbolically from the variational form:

from fenics import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "P", 1)

u = TrialFunction(V)
v = TestFunction(V)

f = Constant(1.0)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx

u_solution = Function(V)
solve(a == L, u_solution)

# Nonlinear: Newton iterate with automatic Jacobian
U = Function(V)  # current solution
u = TrialFunction(V)

F = dot(grad(U), grad(v)) * dx - f * v * dx  # residual
A = derivative(F, U, u)  # Jacobian (symbolic)

problem = NonlinearProblem(F, U, bcs)
solver = NewtonSolver(MPI.comm.world)
solver.parameters["linear_solver"] = "petsc"
solver.solve(problem)

The derivative(F, U, u) call generates the exact Jacobian through automatic differentiation. This is the most powerful feature of FEniCS for nonlinear problems.

PETSc/SNES: Newton-Krylov with Preconditioning

PETSc’s SNES package implements the full Newton-Krylov pattern:

from pysns import SNES

# Define residual function
def residual(x):
    # Return F(x) — the nonlinear residual
    return F_of_x(x)

# Create SNES solver
snes = SNESCreate()
snesSetFunction(snes, x, residual)

# Set Newton-Krylov with preconditioning
SNESSetType(snes, SNESNEWTONKRYLOV)
SNESSetKrylovDimension(snes, 50)  # GMRES max iterations

# Configure preconditioner (physics-based)
ksp = SNESGetKSP(snes)
KSPSetType(ksp, KSPPRECONDEL)  # ILU preconditioner
KSPSetPreconditioner(ksp, PCILU)

# Solve
SNESSetUp(snes)
SNESSolve(snes)

The Newton-Krylov type automatically switches between full Newton (assembled Jacobian), Jacobian-free Newton, and Broyden-based variants depending on the SNESSetUser Jacobian configuration. The preconditioner setup controls the Krylov convergence rate.

NonlinearSolve.jl: The Polyalgorithm Pattern

Julia’s DifferentialEquations ecosystem uses a polyalgorithm that chains methods:

using DifferentialEquations

# Define the nonlinear system
f! = (residual, x) -> begin
    residual[1] = x[1]^2 + x[2] - 1
    residual[2] = x[1] + x[2]^2 - 2
end

# Newton-Krylov with Eisenstein-Walker forcing
prob = NonlinearProblem(f!, [0.5, 0.5])
sol = solve(prob, NewtonRaphson(); abstol=1e-8, reltol=1e-6)

# For large-scale systems: Jacobian-free Newton-Krylov
sol = solve(prob, JFNK(); abstol=1e-8)

The JFNK algorithm computes Jacobian-vector products through finite differences and delegates to a Krylov subspace solver. For stiff systems, the polyalgorithm falls back to TrustRegion if JFNK stalls.

Empirical Performance: JFNK vs Operator-Splitting

A concrete benchmark from the literature demonstrates JFNK’s practical advantage. In non-LTE (NLTE) radiative transfer calculations, a 2024 study in Astronomy & Astrophysics (A&A) compared Jacobian-free Newton-Krylov against operator-splitting approaches. JFNK converged approximately 2× faster than operator-splitting for the same accuracy target, while using comparable memory. The result highlights that monolithic nonlinear solvers can outperform partitioned approaches when the nonlinearity is strong and the preconditioner is well-matched to the physics.

Summary: Which Method Do We Recommend?

The recommendation depends on the scale and structure of your problem:

• Small problems with exact Jacobian available: Newton-Raphson. The quadratic convergence and robust convergence basin outweigh the Jacobian cost.
• Mildly nonlinear 1D-2D systems: Broyden’s method is cost-effective. Accept the risk of condition drift for small-to-medium problems where each iteration is cheap.
• Large-scale 3D systems with millions of degrees of freedom: JFNK with physics-based preconditioning. This is the only viable option when full Jacobian assembly is prohibitive and Broyden’s condition drift is unacceptable.
• Black-box simulation codes: JFNK or Broyden. Both avoid explicit Jacobian construction.
• Critical simulation requirements: Newton with line-search globalization. The quadratic convergence and global stability of trust-region methods are worth the extra cost.

For most production PDE workflows, the polyalgorithm pattern — fast Broyden fallback to Newton fallback to TrustRegion — provides the best balance of speed, robustness, and reliability.

Related Guides

• Time Integration Methods for PDE Solvers: Explicit vs Implicit Schemes — Foundational guide to implicit schemes that require nonlinear solvers.
• Implicit vs Explicit Methods: Stability, Accuracy, and When to Use Each — Discusses stability regions and when implicit methods become necessary.
• Choosing the Right Python PDE Solver: FiPy vs py-pde vs FEniCS — Compares solver libraries, including their nonlinear solver backends.
• Managing Large-Scale PDE Problems: Strategies, Solvers, and HPC Case Studies — Covers nonlinear solvers at scale with HPC infrastructure.
• Operator Splitting Methods: Strang Splitting and IMEX Schemes — Alternatives to monolithic nonlinear solving.

References

• Burden, R., & Faires, J. D. (2011). Numerical Analysis (9th ed.). Cengage Learning. Chapter 8 covers Newton’s method and quasi-Newton methods.
• Wikipedia: Broyden’s method — Overview of the quasi-Newton root-finding family.
• SciPy: scipy.optimize.broyden1 — Documentation for SciPy’s Broyden solver.
• PETSc SNES Documentation — SNU Nonlinear Equation Solver for Newton-Krylov patterns.
• FEniCS Nonlinear Solvers — Automatic Jacobian generation through symbolic differentiation.
• DifferentialEquations.jl: JFNK — Jacobian-free Newton-Krylov for Julia-based PDE solvers.

Want Help Implementing a Nonlinear Solver for Your Research?

Choosing the right nonlinear solver method can make the difference between a simulation that finishes in hours and one that never converges. If you are struggling with Newton iteration divergence, Broyden condition drift, or JFNK preconditioner tuning, our team can help.

We specialize in building verification frameworks for scientific Python and Julia codes, including nonlinear solver workflows for PDE simulations. Get in touch via our issue tracking system to discuss your project’s needs.