calculus and analysis - Numerical partial derivative

For a one-variable numerical function, it's simple to calculate the derivative at a point with Derivative as Szabolcs has pointed out before:

f[x_?NumericQ] := x^2
f'[3.]
(* 6. *)

But this fails for partial derivatives:

g[x_?NumericQ, y_?NumericQ, z_?NumericQ] = x y z + x^2 y^2 z


Derivative[1, 0, 0][g][1., 1., 1.]
(* 3. *)

Derivative[1, 1, 1][g][2., 3., 4.]
(* Unevaluated: Derivative[1, 1, 1][g][2., 3., 4.] *)

ND seems to only handle the one-dimensional case.

Using SeriesCoefficient simply returns the (scaled) Derivative expression:

SeriesCoefficient[g[x, y, z], {x, 2., 1}, {y, 3., 1}, {z, 4., 1}]

(* Derivative[1, 1, 1][g][2., 3., 4.] *)

I'd prefer not to clutter my code with finite difference formulas, since this functionality must be in Mathematica somewhere; where is it?

EDIT: The closest I've found so far is NDSolve`FiniteDifferenceDerivative, but that works on grids and it's a hassle to use for other purposes. Anyone know of a convenient C/Java library that links well with Mathematica and handles all kinds of numerical differentiation?

EDIT2: Does Derivative have accuracy control? (step size or anything)

Clear @ f
f[x_?NumericQ] = Exp[x];
Array[Abs[Derivative[#1][f][1.] - E] &, {8}]

$$\left( \begin{array}{c|c|c} \text{n} & \text{seconds} & \text{error} \\ \hline 5 & 0.123 & 6.77\times 10^{-7} \\ 6 & 0.297 & 0.0000484 \\ 7 & 0.592 & 0.0127 \\ 8 & 1.05 & 1.11 \\ \end{array} \right)$$

Answer

I see no fundamental problem in using ND to answer all your questions. First I'll repeat the definition of your example function, then I do a single and a third partial derivative. Following that, I'll repeat the test of the accuracy for the exponential function:

g[x_?NumericQ, y_?NumericQ, z_?NumericQ] = x y z + x^2 y^2 z

(* ==> x y z + x^2 y^2 z *)

Needs["NumericalCalculus`"]

ND[g[x, 1, 1], x, 1]


(* ==> 3. *)

ND[ND[ND[g[x, y, z], x, 1], y, 1], z, 1]

(* ==> 5. *)

Clear@f
f[x_?NumericQ] = Exp[x];
Array[Abs[
   ND[f[x], {x, #}, 1, WorkingPrecision -> 40, Terms -> 10] - 

    E] &, {8}]

(*
==> {2.29368416218483*10^-21, 9.0878860135398*10^-19, 
 3.069047503987*10^-17, 3.9592354955*10^-16, 3.03377341*10^-15, 
 1.671999*10^-14, 7.334*10^-14, 2.7*10^-13}
*)

The last example with the exponential function is actually discussed specifically in the help for ND, and I just copied the settings from that application.

Edit

With the nested ND calls above, the number of evaluations of the function g may become prohibitively large. Here is a way to reduce the number of derivative evaluations dramatically when doing repeated partial derivatives with ND:

Clear[g, g1, g2];
g[x_?NumericQ, y_?NumericQ, z_?NumericQ] := (c += 1; x y z + x^2 y^2 z)
g1[x_?NumericQ, y_?NumericQ, z_?NumericQ] := ND[g[x1, y, z], x1, x]
g2[x_?NumericQ, y_?NumericQ, z_?NumericQ] := ND[g1[x, y1, z], y1, y]

c = 0;
(* ==> 0 *)

ND[g2[1, 1, z], z, 1] // N

(* ==> 5. *)

c
(* ==> 512 *)

The variable c is just a counter that gets incremented whenever the original function g is called. Compared to ND[Nd[Nd[...]]], the reduction factor is 256.