Skip to main content

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.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more