4. Local sensitivity analysis

Local sensitivity analysis answers a deceptively simple question: if I perturb one design variable by a small amount, how much does a response of interest change? In CAE-driven design optimization, this question is asked thousands of times per iteration, and the way you compute the gradient determines whether an optimization run finishes in minutes or days. This chapter examines the three dominant gradient-computation strategies used inside and around CAE solvers — finite difference, complex step, and adjoint — and explains when each is the right tool.

What a Local Gradient Actually Represents

Consider a response f(x) produced by a solver, where x = [x₁, x₂, …, xₙ] is the vector of design variables (plate thickness, fillet radius, material modulus, and so on). The local sensitivity is the partial derivative ∂f/∂xᵢ evaluated at the current design point. Because the solver is a black box — a sequence of element assembly, linear solve, and post-processing — these derivatives are rarely available in closed form. All three methods below approximate or reconstruct them numerically, but with fundamentally different cost and accuracy profiles.

A critical convention must be fixed before any comparison: throughout this chapter, the perturbation is applied in the positive direction of each design variable, and the gradient sign follows the standard calculus convention (positive means the response increases as the variable increases). Any switch of reference direction will be stated explicitly.

Finite Difference: The Baseline

The forward finite difference approximates the derivative as

∂f/∂xᵢ ≈ [ f(xᵢ + Δ) − f(xᵢ) ] / Δ

and the central difference uses

∂f/∂xᵢ ≈ [ f(xᵢ + Δ) − f(xᵢ − Δ) ] / (2Δ)

The appeal is universality: any solver that accepts an input deck and returns a scalar response can be differentiated this way, with no source-code access. The cost is brutal. Forward differences require n+1 solver runs per gradient; central differences require 2n. For a 40-variable thickness optimization on a model that takes 90 seconds per solve, one central-difference gradient costs 120 solves — roughly three hours — before any optimization logic runs.

Accuracy is the second problem. Finite difference suffers from subtractive cancellation: as Δ shrinks, f(xᵢ+Δ) and f(xᵢ) become nearly equal, and floating-point round-off dominates the numerator. As Δ grows, truncation error from the neglected higher-order terms dominates. The total error therefore has a U-shaped dependence on Δ, with an optimum near Δ ≈ √(ε_mach)·|xᵢ| for forward differences (ε_mach ≈ 2.2×10⁻¹⁶ in double precision). In practice, a step of Δ = 10⁻⁴·|xᵢ| is a reasonable starting point for smooth responses, but noisy responses — contact status changes, mesh-dependent stress peaks, convergence-limited iterative solves — can produce gradients that are essentially garbage at any step size.

Complex Step: Accuracy Without Cancellation

The complex-step method exploits the Taylor expansion of an analytic function along the imaginary axis:

f(x + i·h) = f(x) + i·h·f′(x) − (h²/2)·f″(x) − i·(h³/6)·f‴(x) + …

Taking the imaginary part and dividing by h gives

f′(x) ≈ Im[ f(x + i·h) ] / h

with truncation error O(h²) and — crucially — no subtraction, so no cancellation error. Step sizes as small as h = 10⁻²⁰ can be used, yielding derivatives accurate to machine precision. The catch is that the entire solver must be able to accept complex-valued inputs and propagate them through every operation: element stiffness assembly, the linear solver, material models, and post-processing. Any operation that is not complex-analytic — taking an absolute value, a max/min, a norm, or a real-only branch — silently corrupts the imaginary part and destroys the gradient.

In practice, complex step is implemented by templating the solver's scalar type. A typical C++ pattern looks like:

template <typename T>
T compute_response(const std::vector<T>& x);

// gradient extraction
std::vector<double> grad(n);
for (int i = 0; i < n; ++i) {
    std::vector<std::complex<double>> xc(x.begin(), x.end());
    xc[i] += std::complex<double>(0.0, 1e-20);
    auto fc = compute_response(xc);
    grad[i] = fc.imag() / 1e-20;
}

Cost is still n solves per gradient — the same order as forward differences — but the accuracy is dramatically better, and the method is far more robust to solver convergence noise because the perturbation is infinitesimal. Complex step is the method of choice when you control the solver source and the response is smooth.

Adjoint Method: Cost Independent of Variable Count

For large n, neither finite difference nor complex step scales acceptably. The adjoint method restructures the calculation so that the gradient cost is proportional to the number of responses, not the number of design variables. For a single response and any number of variables, one adjoint solve suffices.

The derivation starts from the discretized equilibrium K(x)·u = F(x), where u is the displacement vector. Differentiating a response f(x, u(x)) gives

df/dxᵢ = ∂f/∂xᵢ + (∂f/∂u)ᵀ · (du/dxᵢ)

The direct approach solves for du/dxᵢ for each variable — again n solves. The adjoint approach instead introduces a Lagrange multiplier λ satisfying

Kᵀ · λ = (∂f/∂u)ᵀ

and then evaluates

df/dxᵢ = ∂f/∂xᵢ − λᵀ · (∂K/∂xᵢ · u − ∂F/∂xᵢ)

The adjoint solve is performed once per response, and the remaining terms are cheap matrix-vector products per variable. For a 40-variable, single-response problem, the adjoint gradient costs one extra linear solve instead of 40 — a 40× reduction. For 500 variables, the advantage is two orders of magnitude.

The price is implementation complexity. The adjoint requires access to the assembled stiffness matrix (or at least its transpose-action), the derivative of the stiffness with respect to each design variable (∂K/∂xᵢ), and the partial derivatives of the response. In commercial solvers this is exposed through dedicated sensitivity or optimization modules; in open-source codes it must be hand-implemented, and the ∂K/∂xᵢ terms must be derived for every element type and design-variable mapping.

Choosing the Method in Practice

The decision is driven by three factors: variable count, solver access, and response smoothness. The table below summarizes the practical trade-offs.

A robust workflow uses more than one method. Compute the gradient with the production method (adjoint or complex step), then verify a few components against central finite differences at a well-chosen step size. If the two disagree by more than a few percent, investigate before feeding the gradient to the optimizer — a bad gradient will stall or misdirect the entire optimization.

Worked Example: Cantilever Tip Deflection

Take a cantilever beam discretized into 20 elements, with element thicknesses as design variables (n = 20) and tip deflection as the single response. A forward finite-difference gradient requires 21 solves; central requires 40; complex step requires 21; adjoint requires 1 adjoint solve plus 20 cheap matrix-vector products. On a solver taking 2 seconds per solve, the adjoint gradient costs roughly 2.5 seconds versus 42 seconds for forward differences — and the gap widens linearly with n. For this problem, the adjoint is clearly the right production method, with central differences reserved for spot-checking two or three components.

The key takeaway is not that one method is universally best, but that the choice must be made deliberately. Finite difference is the universal fallback and the verification standard; complex step is the accuracy champion when the solver can be templated; the adjoint is the scalability champion when variable counts grow. Matching the method to the problem structure is what makes local sensitivity analysis a practical engineering tool rather than a computational bottleneck.