17
$\begingroup$

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

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

2 Answers 2

Reset to default
10
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ 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(!) $\endgroup$ Commented Jun 8, 2013 at 18:29
  • $\begingroup$ Yes, numerical differentiation is known to be sensitive to cancellation errors, so this is to be expected. $\endgroup$ Commented Jun 8, 2013 at 19:01
  • $\begingroup$ @ssch See my edited answer, it's a lot faster... $\endgroup$ Commented Jun 8, 2013 at 22:08
  • $\begingroup$ @Jens, does MMA derivative and ACEGEN derivative work better? $\endgroup$ Commented Mar 19, 2019 at 1:08
  • $\begingroup$ @ABCDEMMM It's certainly an alternative. Maybe it helps to look at the discussion below my answer here where Automatic Differentiation is mentioned. $\endgroup$ Commented Mar 19, 2019 at 3:47
2
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ Nice, sadly I don't have it in symbolic form :( $\endgroup$ Commented Mar 30, 2013 at 21:57
  • $\begingroup$ What form do you have g in? $\endgroup$ Commented Mar 30, 2013 at 23:10
  • $\begingroup$ 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. $\endgroup$ Commented Mar 31, 2013 at 0:28

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.