# Higher-order partial derivatives w.r.t. variable raised to some power

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.

• Maybe read up on Nest? – Marius Ladegård Meyer Jan 19 '17 at 15:30

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 = ReduceFreeVariables[xv];
If[Length[xl] > 1, Return[$Failed]]; fl = ReduceFreeVariables[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.) • 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? – Casimir Jan 21 '17 at 13:18 • 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. – J. M.'s discontentment Jan 21 '17 at 13:23 • Kinda wish I knew what was broken with the method given here so I can fix it; it seems someone was unsatisfied. – J. M.'s discontentment 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  2$\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, ReduceFreeVariables[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]
];
$Failed ) chainD::dvar = "Multiple derivative specifier 1 must be a non-negative integer";  And, here are some examples: chainD[x^4,{x^2,2}] Simplify@chainD[x^5/Sqrt[x^3+y^2],{x^3,2}]//TeXForm  2$\frac{7 x^6+20 x^3 y^2+40 y^4}{36 x \left(x^3+y^2\right)^{5/2}}$Simplify@chainD[(x^3+y^2)/Sqrt[x^3+y^3],x^2,{Log[y],2}]//TeXForm $\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. Basically the same answer I gave you before: Simplify[D[ 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). 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 *)  Then Nest[myD[#, x^2] &, x^4, 2]  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}}$and 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}}$• This works great! – Casimir Jan 20 '17 at 12:22 • 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. – Casimir Jan 20 '17 at 12:23
• 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. – corey979 Jan 20 '17 at 12:56