# 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 NDSolveFiniteDifferenceDerivative, 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)$$

• as an aside, the code for ND is freely visible (it's in AddOns/Packages/NumericalCalculus/NLimit.m).
– acl
Commented Mar 30, 2013 at 17:09
• Related question: How can I differentiate Numerically?
– Jens
Commented May 1, 2013 at 4:07

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.

• The problem with nesting ND is that it becomes incredibly inefficient, ND[ND[ND[g[x, y, z], x, 1], y, 1], z, 1] evaluates the function 131072 times(!)
– ssch
Commented Jun 8, 2013 at 18:29
• Yes, numerical differentiation is known to be sensitive to cancellation errors, so this is to be expected.
– Jens
Commented Jun 8, 2013 at 19:01
• @ssch See my edited answer, it's a lot faster...
– Jens
Commented Jun 8, 2013 at 22:08
• @Jens, does MMA derivative and ACEGEN derivative work better? Commented Mar 19, 2019 at 1:08
• @ABCDEMMM It's certainly an alternative. Maybe it helps to look at the discussion below my answer here where Automatic Differentiation is mentioned.
– Jens
Commented Mar 19, 2019 at 3:47

Here's one way. You have a symbolic base function and numeric top-level one.

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

Derivative[nn__][g][x_?NumericQ, y_?NumericQ, z_?NumericQ] := Derivative[nn][g0][x, y, z]

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

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


One caveat: Somehow the rule is associated to Derivative not g. Clearing g doesn't unset the rule:

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


Clearing Derivative does:

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


(No complaints from Mathematica, and Derivative still works.)

• Nice, sadly I don't have it in symbolic form :(
– ssch
Commented Mar 30, 2013 at 21:57
• What form do you have g in? Commented Mar 30, 2013 at 23:10
• In this case a linked in C function. Sometimes I have symbolic ones where D is prohibitively expensive (minutes) although in those cases I often get by with computing it once and storing that expression to plug numerical values in when needed.
– ssch
Commented Mar 31, 2013 at 0:28