Jacobian-Free Newton-Krylov

Introduction

The Jacobian-Free Newton–Krylov (JFNK) method is a nonlinear solver that combines Newton’s method with Krylov-subspace linear solvers without ever forming the Jacobian matrix explicitly.

At each Newton iteration we solve the linear system

\[J(x^k)\,\delta x = -F(x^k)\]

using a Krylov method (e.g. GMRES). Rather than assembling $J$, JFNK approximates Jacobian–vector products $J(x)v$ via finite differences:

\[J(x)\,v \approx \frac{F(x + \varepsilon\,v) - F(x)}{\varepsilon},\]

with a small $ε$. This “matrix-free” approach reduces memory and computational cost when $J$ is large or expensive to form, and allows reuse of existing residual-evaluation code. Preconditioning can still be applied by approximating $J$ or its inverse action. JFNK is widely used in large-scale simulations (e.g., fluid dynamics, reactive flows) where explicit Jacobian assembly would be prohibitive.

Julia implementation

using LinearAlgebra, LinearMaps, Krylov, Printf

F(x) = [x[1] + x[2]^2, x[1] + 2x[2] + x[1] * x[2]]

function JFNK(F, x0; tol=1e-10, k_max=20, verbose=true)
    n = length(x0)
    x = copy(x0)
    r = copy(F(x))
    # *using `LinearMaps` create a matrix-free operator for Jacobian-vector product J*v ≈ (F(x+δ*v) - F(x)) / δ
    Jop = LinearMap(v -> Jv_product(v, F, x, r), n, n)

    k = 0
    solved = false
    while !solved && k < k_max
        k += 1
        v, stats = Krylov.gmres(Jop, -r)
        x .+= v
        r .= F(x)
        normr = norm(r)

        verbose && @printf(" %2d: x = [%.2e %.2e], ||r|| = %.2e\n", k, x[1], x[2], normr)

        if normr < tol
            if verbose
                printstyled("Equation solved.\n"; color=:green)
                println("Because the residual norm is less than specified tolerance.\n")
            end
            solved = true
        end
        if k == k_max
            if verbose
                printstyled("Failed to solve the equation.\n"; color=9)
                println("JFNK not converge in $k_max iterations.\n")
            end
        end
    end

    return x
end

function Jv_product(v, F, x, r)
    # Finite-difference parameter
    δ = sqrt(eps()) * max(norm(x), 1.0)
    return (F(x .+ δ .* v) .- r) ./ δ
end

x0 = [1.5, 0.1]
x = JFNK(F, x0)

Inputs

• F (Function):
  The nonlinear residual function $F: \mathbb{R}^n \to \mathbb{R}^n$. Each call F(x) returns an $n$-element vector.

• x0 (AbstractVector{<:Real}):
  The initial guess for the solution, of length $n$. The algorithm will refine this vector.

• tol (Real, keyword; default 1e-10):
  The stopping tolerance on the residual norm. Once $|F(x)| < \text{tol}$, iteration stops.

• k_max (Integer, keyword; default 20):
  The maximum number of Newton–Krylov (outer) iterations allowed before giving up.

• verbose (Bool, keyword; default true):
  If true, print progress each iteration (values of x and $|r|$) and a final success/failure message with colored output.

Internal Variables

• n: problem size (length of x0).

• x: current iterate (updated each iteration).

• r: current residual, $r = F(x)$.