Applying Gradient Operator on Nonuniform Grids

What Is a Gradient Operator and Why Grid Type Matters

The gradient operator calculates spatial derivatives. In physics and engineering, this means finding how a quantity changes across space—temperature gradients, pressure fields, velocity profiles. Pretty much every simulation you run depends on getting this right.

Most textbooks show gradient calculations on uniform grids where spacing between points is constant. That's clean. That's easy. That's also rarely what you have in real applications.

Nonuniform grids appear constantly:

• Adaptive mesh refinement where you need finer resolution near boundaries
• Geophysical grids following Earth's curvature
• Computational fluid dynamics with boundary layer clustering
• Any problem with natural scaling—financial models, population distributions

When your grid points aren't evenly spaced, standard finite difference formulas break down. You need different math.

The Core Problem with Nonuniform Spacing

On a uniform grid, the first derivative approximates as:

f'(x) ≈ (f(x+h) - f(x-h)) / (2h)

This assumes h is constant everywhere. When spacing varies, that assumption fails. The error terms change, and your approximation degrades or becomes outright wrong.

The gradient operator on nonuniform grids requires accounting for the actual spacing between each pair of points. This isn't optional tweaking—it's fundamental to getting correct results.

Methods for Computing Gradients on Nonuniform Grids

Finite Difference Methods

The most common approach. You derive coefficients based on actual point positions rather than assuming uniform spacing.

For two adjacent points xᵢ and xᵢ₊₁ with spacing hᵢ = xᵢ₊₁ - xᵢ, a first-order backward difference gives:

f'(xᵢ) ≈ (fᵢ - fᵢ₋₁) / hᵢ₋₁

Second-order accurate central differences on nonuniform grids use:

f'(xᵢ) ≈ [hᵢ² fᵢ₋₁ - (hᵢ₋₁² - hᵢ²) fᵢ - hᵢ₋₁² fᵢ₊₁] / [hᵢ₋₁ hᵢ (hᵢ₋₁ - hᵢ)]

This looks ugly, but it's just algebra. The coefficients depend on your local spacing, nothing more.

Polynomial Interpolation Approach

Fit a polynomial through nearby grid points, then differentiate analytically. More flexible but computationally heavier.

Lagrange interpolation is the standard choice. For three points, the derivative at xᵢ becomes:

f'(xᵢ) ≈ Σⱼ fⱼ · Lⱼ'(xᵢ)

where Lⱼ are Lagrange basis polynomials evaluated at your grid points. The derivative of each basis polynomial accounts for the nonuniform spacing automatically.

Spectral Methods

For smooth solutions on structured nonuniform grids, spectral collocation works. You map your physical grid to a uniform computational space, apply spectral differentiation there, then transform back.

This works well when your grid is smoothly varying. It falls apart when you have sharp transitions or discontinuities.

Comparison of Methods

Getting Started: Practical Implementation

Here's how to compute gradient operators on nonuniform grids in Python. No black-box libraries—just the math.

Second-Order Central Difference on Nonuniform Grid

import numpy as np

def gradient_nonuniform(f, x):
    """
    Compute gradient using 2nd-order central difference.
    Works on arbitrary nonuniform grids.
    
    Parameters:
        f: array of function values at grid points
        x: array of grid point positions
    
    Returns:
        dfdx: gradient at interior points
    """
    n = len(x)
    dfdx = np.zeros(n)
    
    for i in range(1, n-1):
        h_left = x[i] - x[i-1]
        h_right = x[i+1] - x[i]
        
        # Second-order accurate central difference
        dfdx[i] = (h_right**2 * f[i-1] 
                   - (h_left**2 - h_right**2) * f[i] 
                   - h_left**2 * f[i+1]) / (h_left * h_right * (h_left - h_right))
    
    return dfdx

Handling Boundaries

Central differences don't work at boundaries. Use one-sided differences there:

def gradient_with_boundaries(f, x):
    """
    Full gradient computation with proper boundary handling.
    """
    n = len(x)
    dfdx = np.zeros(n)
    
    # Interior points: central difference
    dfdx[1:-1] = gradient_nonuniform(f[1:-1], x[1:-1])
    
    # Left boundary: forward difference
    h0 = x[1] - x[0]
    dfdx[0] = (-f[2] + 4*f[1] - 3*f[0]) / (2*h0)
    
    # Right boundary: backward difference  
    h_last = x[-1] - x[-2]
    dfdx[-1] = (f[-3] - 4*f[-2] + 3*f[-1]) / (2*h_last)
    
    return dfdx

Testing Your Implementation

Always verify against known analytical results. Here's a test case:

# Test with known function
x = np.array([0, 0.3, 0.7, 1.2, 1.8, 2.5])  # Nonuniform
f = x**3  # f(x) = x^3, f'(x) = 3x^2

dfdx_computed = gradient_with_boundaries(f, x)
dfdx_exact = 3 * x**2

print("Computed:", dfdx_computed)
print("Exact:   ", dfdx_exact)
print("Error:   ", np.abs(dfdx_computed - dfdx_exact))

For smooth functions, expect errors in the 10⁻⁶ to 10⁻⁸ range with double precision. Larger errors mean your implementation has bugs or you're using insufficient order of accuracy.

Common Pitfalls

• Ignoring boundary treatment — Central schemes don't apply at edges. Forgetting this introduces large errors that contaminate your entire solution.
• Assuming smoothly varying spacing — Some methods break down when grid spacing changes abruptly. Check if your method handles this.
• Mixing coordinate systems — In curvilinear coordinates, the gradient includes metric terms. The formulas above assume Cartesian coordinates.
• Machine precision issues — When grid points are extremely close together, roundoff errors dominate. Watch your condition numbers.

When to Use Which Method

For most engineering problems, second-order finite differences on nonuniform grids give the best tradeoff between accuracy and simplicity. They're easy to implement, fast to compute, and work on essentially any grid topology.

Use Lagrange interpolation when you need higher accuracy and have smooth solutions. Increase the stencil size to gain accuracy at the cost of more computations.

Use spectral methods only when you have smoothly varying grids and genuinely smooth solutions. The moment you have discontinuities or sharp gradients, spectral methods produce Gibbs oscillations that ruin your results.

The gradient operator on nonuniform grids isn't complicated once you see that the formulas are just accounting for local spacing. The math isn't harder—it's the same derivatives, just with spacing variables instead of constants. Implement once, test thoroughly, and reuse everywhere.