Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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)$$

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

2 Answers

up vote 4 down vote accepted

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.

share|improve this answer
    
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 Jun 8 '13 at 18:29
    
Yes, numerical differentiation is known to be sensitive to cancellation errors, so this is to be expected. –  Jens Jun 8 '13 at 19:01
    
@ssch See my edited answer, it's a lot faster... –  Jens Jun 8 '13 at 22:08
add comment

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.)

share|improve this answer
    
Nice, sadly I don't have it in symbolic form :( –  ssch Mar 30 '13 at 21:57
    
What form do you have g in? –  Michael E2 Mar 30 '13 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 Mar 31 '13 at 0:28
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.