# Recursive Formal Differentiation with Substitution

I have a function f[t] and I'm really hoping this question has an answer treating this simply as a formal function, and not requiring the input of its closed form. We know the function satisfies the following two constraints

(D[f[t],t])^2 = 4*(f[t])^3 - A*f[t] - B;
D[f[t],{t,2}] = 6*(f[t])^2 - (1/2)*A;


Here, A and B are constants with respect to t.

(Some may have now noticed that I'm talking about the Weierstrass $\wp$-function $\wp(\tau, z)$ where the derivatives are with respect to $z$. However, perhaps I'm wrong, but I don't think I want to use the built-in Mathematica function for $\wp$. I would like to be more formal)

Now, I'm interested in the even number of derivatives D[f[t],{t,2n}]. It's easy to convince oneself that such an even number of derivatives can be expressed as a polynomial in f, A, and B. The idea is that starting from the second of the two equations above, you can differentiate with respect to t an even number of times and you can substitute both of the above equations to get rid of the first derivative squared and second derivatives.

Quite simply, my question is the following: can someone help me write a code as a function of $n$ which will spit out the polynomial of f, A, and B which equals D[f[t],{t,2n}]? I do not want f[t] to be part of the input; I'd like it to be totally formal or symbolic. This should be able to be done recursively. As I mentioned above, if one is given the $(2n-2)$-th derivative, compute the $2n$-th derivative and then substitute in those two equations I provide where appropriate.

(I'm sorry if this is a hopeless question, I convinced myself it should be possible if I was better at Mathematica!)

This returns the polynomials for the derivatives $0,2,4,\dotsc,d$ all at once ($d = 12$ in the example):

d = 12;
repo1 = (D[f[t], t])^2 -> 4*(f[t])^3 - A*f[t] - B;
repo2 = D[f[t], {t, 2}] -> 6*(f[t])^2 - (1/2)*A;
NestList[Expand[D[#, {t, 2}] /. repo2 /. repo1] &, f[t], Quotient[d,2]] /. f[t] -> f

• Please, revise your code. Even correcting your typos, a check for $n=4$ it seem it does not agree with a manual calculation... – José Antonio Díaz Navas Apr 27 '18 at 11:04
• Thanks for the edits. That haven't been typos; OP changed the formlas after I answered. I added also some short description of the return value. Maybe that resolves the other issue? – Henrik Schumacher Apr 27 '18 at 11:08
• OK, now agrees with my proposal ;))... – José Antonio Díaz Navas Apr 27 '18 at 11:14
• I got your results for $d=14$ ??? and missed $d=0$ which should be $f(t)$ :(( ... – José Antonio Díaz Navas Apr 27 '18 at 11:17
• Ah, thank you. Very good point! Fixed. – Henrik Schumacher Apr 27 '18 at 11:21

This also provides the requested result:

r2 = Derivative[2][f][t] -> 6*(f[t])^2 - (1/2)*A;
r1 = Power[(Derivative[1][f][t]), n_] -> Power[Sqrt[4*(f[t])^3 - A*f[t] - B], n];
dn[n_] := Expand@Last@NestList[(D[#, {t, 1}] //. {r1, r2}) &, f[t], n] /. r1

(Expand@dn[#] & /@ Range[2, 8, 2]) // TableForm


$6 f(t)^2-\frac{A}{2}\\-18 A f(t)-12 B+120 f(t)^3\\9 A^2-1008 A f(t)^2-720 B f(t)+5040 f(t)^4\\3024 A^2 f(t)+2376 A B-90720 A f(t)^3-64800 B f(t)^2+362880 f(t)^5$

• Thanks a lot! Although this gave an error when I first ran it. I believe you mean to have an r2 after rule2? That seemed to clear it up on my end. – Benighted Apr 27 '18 at 16:37
• Oh, yes, well spotted !! Corrected... – José Antonio Díaz Navas Apr 27 '18 at 16:38
• Nice, thanks a lot. I'd accept both answers if I could! I appreciate your time. – Benighted Apr 27 '18 at 20:45

Here's an alternate approach. Define DownValues for derivatives of f (and of an auxiliary function g):

Derivative[n_?Positive][f] ^:= fExpand @ Derivative[1][Derivative[n-1][f]]
Derivative[1][f] ^= g;

Derivative[n_?Positive][g] ^:= fExpand @ Derivative[1][Derivative[n-1][g]]
Derivative[1][g] ^= Function[t, 6f[t]^2-A/2];

Power[g[t], n_?(GreaterThan[1])] ^:= (4f[t]^3-A f[t]-B)^Quotient[n,2] g[t]^Mod[n,2]


The auxiliary function g satisfies $g(t)^2 = 4 f(t)^3 - A f(t) - B$ so that the first constraint is satisfied, while the fExpand function is used to reduce intermediate expression swell. Let's check:

f''[t]
Derivative[4][f][t]
Derivative[6][f][t]
Derivative[8][f][t]


-(A/2) + 6 f[t]^2

-12 B - 18 A f[t] + 120 f[t]^3

9 A^2 - 720 B f[t] - 1008 A f[t]^2 + 5040 f[t]^4

2376 A B + 3024 A^2 f[t] - 64800 B f[t]^2 - 90720 A f[t]^3 + 362880 f[t]^5

Derivative[22][f][t] //AbsoluteTiming