In the previous post, Complex Magic, Part 1, the use of complex variables in approximating the *first* derivative of a holomorphic function was discussed. This post delves deeper and introduces formulas suitable for approximating the th derivative of a holomorphic function.

As elaborate as they might seem, the ideas to be discussed are based on two familiar topics: Cauchy’s *integral* theorem and quadratures.

Cauchy’s integral theorem states that if is a holomorphic function in a neighborhood which contains a closed disk with the boundary (a closed rectifiable curve with winding number one about ), then for every point in the interior of that disk

The magic is that (1) converts the problem of calculating a derivative to calculating an integral! And contour integrals of holomorphic functions can be readily approximated using trapezoidal rule, with exponential accuracy. What follows is how this approximation is done, based on [LynessMoler1967].

On the interval partitioned into segments of equal length, the trapezoidal quadrature may be defined with the following sum operator:

Suppose is a real function of whose *n*th derivative we wish to approximate at . If satisfies the conditions required for (1) to hold at , then the derivative can be calculated according to

in which and is the Möbius function, defined for all positive integers as , if is squarefree and even, , if is squarefree and odd, and otherwise. With the sum going to , (3) is, of course, an *exact* formula. Truncating the infinite sum to an -term partial sum yields an *approximate* formula for the derivative.

A condensed version of the proof of [LynessMoler1967] for (3) is as follows: substituting the parametrization in (1) gives the Fourier coefficients

which, using the Poisson summation formula, can be connected to the trapezoidal rule sum operator (2) to yield

The set of equations is then inverted, using the Möbius inversion formula, to yield in which substituting for from (1), with , and for (since for real, is real) completes the proof. Notice that only in the final step of the proof was the fact that is a real function used.

The parametrization , is in effect taking the contour to be a circle of radius centered at . The important point about this parametrization is *not* so much the particular shape that it gives to , but the symmetry of that shape with respect to the real axis. Furthermore, for a real function , is also symmetric with respect to , or equivalently is also symmetric with respect to the real axis. This symmetry of and , together with the circular shape of the former, can be exploited to simplify (2). Specifically, since , (2) reduces simply to a mean of over , which can be approximated by taking the real part of the means of over equally spaced points along the upper (or lower) half of . It is emphasized that because is a real function, the mean of and over (the *full* circle) are equal since the imaginary parts of from the upper and lower halves of cancel each other (up to roundoff).

It was mentioned in Complex Magic, Part 1 that finite difference formulas have coefficients that vary significantly in magnitude, and that this causes the absolute errors in function values to be amplified. Notice that all coefficients in (3) are of magnitude and as a result they will not cause any amplification of the absolute errors in function values. Consequently the partial sum is expected to have an absolute error of the same order as that of the largest function value contributing to it, which is and since dividing by retains the relative accuracy, can serve as a measure of the relative error incurred in calculating the derivative due to roundoff. ( is the maximum relative roundoff error in function evaluations.) If is not acceptable, we can increase the radius, , of the contour , to have the function evaluations done on points prominently distinguishable from (in floating point sense), and repeat the approximation.

Here is a simple implementation of the Möbius function, using table look-up for the first 77 values and factorization for the rest:

function mu = mobius(n) persistent muList if isempty(muList) muList = [1,-1,-1,0,-1,1,-1,0,0,1,-1,0,-1,1,1,0,-1,0,-1,0,... 1,1,-1,0,0,1,0,0,-1,-1,-1,0,1,1,1,0,-1,1,1,0,-1,... -1,-1,0,0,1,-1,0,0,0,1,0,-1,0,1,0,1,1,-1,0,-1,1,0,... 0,1,-1,-1,0,1,-1,-1,0,-1,1,0,0,1]; end if n <= length(muList), mu = muList(n); return, end fac = factor(n); len_fac = length(fac); if length(unique(fac)) == len_fac mu = (-1)^len_fac; else mu = 0; end

And here is a function `diffc`

that implements (3):

function [S_M,delta_M] = diffc(f,varargin) [x0,n,r,M] = deal(0,1,.1,3); optArgs = {x0,n,r,M}; emptyArgs = cellfun(@isempty,varargin); [optArgs{~emptyArgs}] = varargin{~emptyArgs}; [x0,n,r,M] = optArgs{:}; f_x0 = f(x0); [G_M,S_M] = deal(0); for m = 1:M t = linspace(1/(m*n),1,m*n); g = real( f( x0 + r * exp(1i * 2*pi * t) ) ); G_M = max( abs( [G_M f_x0 g] ) ); S_M = S_M + mobius(m) * ( mean(g) - f_x0 ); end delta_M = eps(G_M) * G_M/S_M; S_M = factorial(n)/(r^n) * S_M;

The following code snippet puts `diffc`

to work, approximating the 10th derivative of the function at , against an exact calculation using symbolic variables. The exact value of this derivative is 13829824, as seen from the symbolic calculation.

syms x real f_sym = exp(x)/(sin(x)^3 + cos(x)^3) f = matlabFunction(f_sym) x_0 = 0; n = 10; r = .5; M = 7; tic, df10d10x_sym = subs(diff(f_sym,n),x_0), toc tic, [df10d10x_M,delta_M] = diffc(f,x_0,n,r,M), toc

On a contour of radius using a partial sum of only terms, `diffc`

achieves 13829824.00000018 and the time it takes is two orders of magnitude less than that of the symbolic calculation (0.001s compared to 0.6s).

**Note**

There exists a bug in MATLAB R2009a (5.2) and R2008b (5.1) due to which the “MuPAD kernel can take several minutes or longer to initialize on Windows. Any command from the Symbolic Math Toolbox that requires starting MuPAD might seem to hang.” This bug is fixed in R2009b and a workaround is available for the two releases before that, http://www.mathworks.com/support/bugreports/523801.(When I was first writing this post, without the “patch” applied, I noticed that the symbolic calculation was taking an unreasonably long time, somewhere around 15s. As satisfied as I was as with my uber-fast “complex” implementation, I knew that something was wrong with the Symbolic Math Toolbox because even creating a symbolic variable, say `syms x real`

, would take over 10s.)