Numerical differentiation with a complex step

The most obvious way to approximate the derivative of a function numerically is to use the definition of derivative and stick in a small value of the step size h.

f′ (x) ≈ ( f(x + h) − f(x) ) / h.

How small should h be? Since the exact value of the derivative is the limit as h goes to zero, the smaller h is the better. Except that’s not true in computer arithmetic. When h is too small, f(x + h) is so close to f(x) that the subtraction loses precision.

One way to see this is to think of the extreme case of h so small that f(x + h) equals f(x) to machine precision. Then the derivative is approximated by 0, regardless of what the actual derivative is.

As h gets smaller, the approximation error decreases, but the numerical error increases, and so there is some optimal value of h.

You can do significantly better by using a symmetric approximation for the derivative:

f′ (x) ≈ ( f(x + h) − f(x − h) ) / 2h.

For a given value of h, this formula has about the same numerical error but less approximation error. You still have a trade-off between approximation error and numerical error, but it’s a better trade-off.

If the function f that you want to differentiate is analytic, i.e. differentiable as a function of a complex variable, you can take the step h to be a complex number. When you do, the numerical difficulty completely goes away, and you can take h much smaller.

Suppose h is a small real number and you take

f′ (x) ≈ ( f(x + ih) − f(x − ih) ) / 2ih.

Now f(x + ih) and −f(x − ih) are approximately equal by the Schwarz reflection principle. And so rather than canceling out, the terms in the numerator add. We have

f(x + ih) − f(x − ih) ≈ 2 f(x + ih)

and so

f′ (x) ≈ f(x + ih) / ih.

Example

Let’s take the derivative of the gamma function Γ(x) at 1 using each of the three methods above. The exact value is −γ where γ is the Euler-Mascheroni constant. The following Python code shows the accuracy of each approach.

    from scipy.special import gamma

    def diff1(f, x, h):
        return (f(x + h) - f(x))/h

    def diff2(f, x, h):
        return (f(x + h) - f(x - h))/(2*h)

    γ = 0.5772156649015328606065

    x     = 1
    h     = 1e-7
    exact = -γ

    approx1 = diff1(gamma, x, h)
    approx2 = diff2(gamma, x, h)
    approx3 = diff2(gamma, x, h*1j)

    print(approx1 - exact)
    print(approx2 - exact)
    print(approx3 - exact)

This prints

    9.95565755390615e-08
    1.9161483510998778e-10
    (9.103828801926284e-15-0j)

In this example the symmetric finite difference approximation is about 5 times more accurate than the asymmetric approximation, but the complex step approximation is 10,000,000 times more accurate.

It’s a little awkward that the complex step approximation returns a complex number, albeit one with zero complex part. We can eliminate this problem by using the approximation

f′ (x) ≈ Im f(x + ih) / h

which will return a real number. When we implement this in Python as

    def diff3(f, x, h):
        return (f(x + h*1j) / h).imag

we see that it produces the same result as diff2 but without the zero imaginary part.

Related posts

• Applied complex analysis
• Joukowski irfoils
• Bounding complex roots by real roots