This is a follow-up question to Differentiate w.r.t. squared variable.

I want to be able to compute partial derivatives of the form $$\partial_{x^n}^m \, f(x) = \frac{\partial^m}{\partial^m x^n} \, f(x).$$ For instance, $\partial_{x^2}^2 \, x^4$ should give $2$, and $$\partial_{x^3}^2 \frac{x^5}{\sqrt{x^3+y^2}} = \frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}.$$

What I have so far is

myD[f_, x_^n_] := FullSimplify[D[f, x]/D[x^n, x]]

This works great for computing first-order partial derivatives with respect to quantities raised to some power, i.e. derivatives of the form $\partial_{x^n} \, f(x)$.

To compute higher-order derivatives $\partial_{x^n}^m \, f(x)$ for $m > 1$, I would somehow need to store intermediate results after each differentiation and divide by the factor $\partial x^n/\partial x$ from the chain rule, before applying the next partial derivative. I'm fairly new to Mathematica, and I have no idea how to do this. Any help would be great.


Using a formula adapted from Todorov's paper,

powerD[f_, x_^(k_.)] := powerD[f, {x^k, 1}];
powerD[f_, {x_^(k_.), 0}] := f;
powerD[f_, vars__] := Fold[powerD, f, {vars}];
powerD[f_, {x_^(k_.), n_Integer?Positive}] := 
       Det[Append[Table[(j!/i!) Binomial[k i, j] x^(k i - j), {i, n - 1}, {j, n}], 
                  Table[D[f, {x, j}], {j, n}]]]/(k x^(k - 1))^Binomial[n + 1, 2];

Some examples:

powerD[x^5/Sqrt[x^3 + y^2], x^2] // Simplify
   (x^3 (7 x^3 + 10 y^2))/(4 (x^3 + y^2)^(3/2))

powerD[x^5/Sqrt[x^3 + y^2], {x^3, 2}] // Simplify
   (7 x^6 + 20 x^3 y^2 + 40 y^4)/(36 x (x^3 + y^2)^(5/2))

powerD[x^5/Sqrt[x^3 + y^2], {y^3, 3}] // Simplify
   -((x^5 (4 x^6 + 17 x^3 y^2 + 28 y^4))/(27 y^7 (x^3 + y^2)^(7/2)))

powerD[x^5/Sqrt[x^3 + y^2], {x^2, 2}, {y, 3}] // Simplify
   (45 x (9 x^9 y - 88 x^6 y^3 + 84 x^3 y^5 - 8 y^7))/(16 (x^3 + y^2)^(11/2))

One can even use Todorov's results to define a more general function parametricD[]:

parametricD[f_, xv_] := parametricD[f, {xv, 1}];
parametricD[f_, {xv_, 0}] := f;
parametricD[f_, rest__] := Fold[parametricD, f, {rest}];

parametricD[f_, {xv_, n_Integer?Positive}] := Module[{fl, xl, xx},
          xl = Reduce`FreeVariables[xv];
          If[Length[xl] > 1, Return[$Failed]];
          fl = Reduce`FreeVariables[f];
          If[Intersection[xl, fl] === {}, Return[0]];
          xx = First[xl];
          Det[Append[Table[D[xv^i, {xx, j}]/i!, {i, n - 1}, {j, n}], 
                     Table[D[f, {xx, j}], {j, n}]]]/D[xv, xx]^Binomial[n + 1, 2]]

Some examples:

parametricD[x^5/Sqrt[x^3 + y^2], {x^3, 2}] // Simplify
   (7 x^6 + 20 x^3 y^2 + 40 y^4)/(36 x (x^3 + y^2)^(5/2))

parametricD[(x^3 + y^2)/Sqrt[x^3 + y^3], x^2, {Log[y], 2}] // Simplify
   (3 x y^2 (x^3 (94 - 99 y) y^3 + 2 x^6 (-8 + 9 y) + y^6 (-25 + 18 y)))/
   (16 (x^3 + y^3)^(7/2))

parametricD[x y, BesselJ[2, x], BesselJ[3, y]] // Simplify
   4/((BesselJ[1, x] - BesselJ[3, x]) (BesselJ[2, y] - BesselJ[4, y]))

(Thanks to Mike Spivey for pointing out Todorov's paper.)

|improve this answer|||||
  • $\begingroup$ Works like a charm, but I don't understand how. Could you explain what the last line in your definition of powerD does? For instance, why is Append needed? Also, am I understanding the documentation correctly in that the pattern k_. adopts the default value k = 1 from the first line? $\endgroup$ – Casimir Jan 21 '17 at 13:18
  • 1
    $\begingroup$ I linked to a paper in the very first line; you'll want to read it. The formula came from there. As for k_.: have a look at what Power[x] and DefaultValues[Power] returns. $\endgroup$ – J. M.'s technical difficulties Jan 21 '17 at 13:23
  • $\begingroup$ Kinda wish I knew what was broken with the method given here so I can fix it; it seems someone was unsatisfied. $\endgroup$ – J. M.'s technical difficulties Apr 10 '18 at 3:08

Another method is to use UpValues as follows:

chainD[f_, z:{g_,_}|g_, x_] := Module[{v},
    v /: D[x, v, NonConstants->{x}] = 1/D[g,x];
    D[f, z /. g->v, NonConstants->{x}]

Then we find:

chainD[x^4, {x^2,2}, x]
Simplify @ chainD[x^5/Sqrt[x^3+y^2], {x^3, 2}, x] //TeXForm


$\frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}$

One can refine the method so that the variable x need not be specified, (edit) and to handle multiple derivatives:

chainD[f_, z_, r__] := Module[{res = Check[chainD[f, z], $Failed]},
    chainD[res, r] /; res=!=$Failed

chainD[f_, z_] := Module[{res = iD[f, z, Reduce`FreeVariables[z /. {g_, _} :> g]]},
    res /; res =!= $Failed

iD[f_, {g_, n_Integer?Positive}, {x_}] := Module[{v},
    v /: D[x, v, NonConstants->{x}] = 1/D[g, x];
    D[f, {v, n}, NonConstants->{x}]

iD[f_, g:Except[_List], {x_}] := Module[{v},
    v /: D[x, v, NonConstants->{x}] = 1/D[g,x];
    D[f, v, NonConstants->{x}]

iD[_, z:{g_, n_} | g_, x_] := (
    If[!MatchQ[x, {_}],
        Message[General::ivar, g],
        Message[chainD::dvar, n]

chainD::dvar = "Multiple derivative specifier `1` must be a non-negative integer";

And, here are some examples:



$\frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}$


$\frac{3 x \left(2 x^6 (9 y-8) y^2+x^3 (94-99 y) y^5+(18 y-25) y^8\right)}{16 \left(x^3+y^3\right)^{7/2}}$

Simplify @ chainD[x y, BesselJ[2,x], BesselJ[3,y]]//TeXForm

$\frac{4}{(J_1(x)-J_3(x)) (J_2(y)-J_4(y))}$

Note that contrary to some other suggested answers, chainD works with any function, not just powers or invertible functions.

|improve this answer|||||

Basically the same answer I gave you before:

   x^5/Sqrt[x^3 + y^2] /. x -> asdfljkhlh^(1/3), {asdfljkhlh, 2}] /. 
  asdfljkhlh -> x^3, Assumptions -> x > 0]
(* (7 x^6 + 20 x^3 y^2 + 40 y^4)/(36 x (x^3 + y^2)^(5/2)) *)

The advantage of this concept of replacing and back-replacing is that you can easily calculate higher order derivatives. The naming of the temporary variable is such that you won't run into name conflicts with existing variables (at least it would be highly unlikely that you name your other variables the same).

|improve this answer|||||

Basically the same answer as chuy gave before, with a proper use of Nest:

myD[f_, x_, cons : OptionsPattern[Dt]] := FullSimplify[Dt[f, cons]/Dt[x, cons]]
 (* chuy's function *)


Nest[myD[#, x^2] &, x^4, 2]


Nest[myD[#, x^3, Constants -> y] &, x^5/Sqrt[x^3 + y^2], 2]

$\frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}$

Wrapping in a function:

myDn[func_, var_, n_, cons : OptionsPattern[Dt]] := Nest[myD[#, var, cons] &, func, n]

so that

myDn[x^5/Sqrt[x^3 + y^2], x^3, 2, Constants -> y]

is obviously also

$\frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}$


myDn[x^5/Sqrt[x^3 + y^2], y^3, 3, Constants -> x]

$-\frac{x^5 \left(4 x^6+17 x^3 y^2+28 y^4\right)}{27 y^7 \left(x^3+y^2\right)^{7/2}}$

|improve this answer|||||
  • $\begingroup$ This works great! $\endgroup$ – Casimir Jan 20 '17 at 12:22
  • $\begingroup$ You take derivatives w.r.t. both $x$ and $y$. Is it possible to modify myDn so as to be able to differentiate w.r.t. to multiple variables in one go, i.e. write something like myDn[f,{x^m,n},{y^p,q}]? That would bring myDn up to par with Mathematica's built-in derivative functions. $\endgroup$ – Casimir Jan 20 '17 at 12:23
  • $\begingroup$ While myDn[myDn[x^5/Sqrt[x^3 + y^2], x^2, 2, Constants -> y], y, 3, Constants -> x] works, making it a function is not that straightforward. J.M. has provided a general answer. $\endgroup$ – corey979 Jan 20 '17 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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