Skip to main content

Jacobian-Free Newton-Krylov

· 3 min read

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(xk) δx=−F(xk)J(x^k)\,\delta x = -F(x^k)

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

J(x) v≈F(x+ε v)−F(x)ε,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 JJ is large or expensive to form, and allows reuse of existing residual-evaluation code. Preconditioning can still be applied by approximating JJ 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:Rn→RnF: \mathbb{R}^n \to \mathbb{R}^n. Each call F(x) returns an nn-element vector.

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

  • tol (Real, keyword; default 1e-10): The stopping tolerance on the residual norm. Once ∥F(x)∥<tol\|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∥\|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)r = F(x).