higher order derivatives of Bessel functions

[rtoy]

Imported from SourceForge on 2024-07-09 11:14:58 Created by willisbl on 2023-11-07 16:12:59 Original: https://sourceforge.net/p/maxima/bugs/4207

---

The 25th derivative of bessel_j(0,x) isn't particularly lengthy, but Maxima can't do it:

(%i6)   diff(bessel_j(0,x),x,25)$
Message from the stdout of Maxima: Welcome to LDB,...

This is Maxima 5.47 compiled with SBCL for Windows.

[rtoy]

Imported from SourceForge on 2024-07-09 11:14:59 Created by robert_dodier on 2023-11-07 18:16:04 Original: https://sourceforge.net/p/maxima/bugs/4207/#0778

---

The following seems to show the size of the expression approximately doubles for each additional derivative ...

expr_size (e) := if atom(e) then 1 else 1 + lsum (expr_size (a), a, args (e));
makelist (expr_size (diff (bessel_j (0, x), x, k)), k, 1, 12);
 => [4, 11, 19, 40, 70, 140, 250, 495, 901, 1783, 3295, 6529]

Given that, it's not surprising that sooner or later Maxima runs out of steam.

ratsimp reduces those expressions pretty well; should diff call ratsimp? I don't know what the implications of that would be. Perhaps it's worth a try to investigate at least.

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:02 Created by willisbl on 2023-11-07 18:22:44 Original: https://sourceforge.net/p/maxima/bugs/4207/#518f

---

Agreed--almost surely inserting a call to sratsimp would be the right thing to do

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:06 Created by macrakis on 2023-11-08 08:33:36 Original: https://sourceforge.net/p/maxima/bugs/4207/#518f/5001

---

I disagree. In some cases, ratsimp does make the expression smaller. In other cases, it makes it bigger, e.g.:

(%e^x+x)^5$

diff(%,x);
  5*(%e^x+1)*(%e^x+x)^4

ratsimp(%)
  5*%e^(5*x)+(20*x+5)*%e^(4*x)+(30*x^2+20*x)*%e^(3*x)
                 +(20*x^3+30*x^2)*%e^(2*x)+(5*x^4+20*x^3)*%e^x+5*x^4

Not only bigger, but it has lost the structure.

The user can apply ratsimp or other transformations (expand, factor, etc.) when they find it useful.

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:10 Created by willisbl on 2023-11-08 11:46:54 Original: https://sourceforge.net/p/maxima/bugs/4207/#9013

---

Oh, I see. I was thinking there would be a way to apply ratsimp non-globally (only to the result of each derivative of a Bessel function). I doubt that's possible.

As I recall, long ago, Maxima had a specialized nth derivative function. It was so buggy that it was expunged.

Even if we could squeeze in a call toratsimp or to rat, this code has terrible big-O behavior:

(%i27)  dn(e,x,n) := if n=0 then rat(e) else(rat(dn(diff(e,x),x,n-1)))$

(%i33)  dn(bessel_j(0,x),x,18)$

(%i34)  time(%);
(%o34)  [9.671875]

(%i35)  dn(bessel_j(0,x),x,19)$

(%i36)  time(%);
(%o36)  [20.546875]

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:13 Created by macrakis on 2023-11-08 12:24:04 Original: https://sourceforge.net/p/maxima/bugs/4207/#baa9

---

Surprisingly, perhaps, one big rat/ratsimp at the end costs about the same as incremental rats or ratsimps (final expand is a bit slower):

Rat

dn(e,x,n) := if n=0 then rat(e) else(rat(dn(diff(e,x),x,n-1)))$
dn(bessel_j(0,x),x,18)$
2.8 sec

rat(diff(bessel_j(0,x),x,18))$
2.8 sec

Ratsimp

dn(e,x,n) := if n=0 then rat(e) else(ratsimp(dn(diff(e,x),x,n-1)))$
dn(bessel_j(0,x),x,18)$
2.6 sec

ratsimp(diff(bessel_j(0,x),x,18))$
2.7 sec

Expand

dn(e,x,n) := if n=0 then rat(e) else(expand(dn(diff(e,x),x,n-1)))$
dn(bessel_j(0,x),x,18)$
2.7 sec

expand(diff(bessel_j(0,x),x,18))$
3.7 sec

Maxima 5.46.0 SBCL 2.3.0 MacOS

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:16 Created by macrakis on 2023-11-08 12:44:14 Original: https://sourceforge.net/p/maxima/bugs/4207/#c15f

---

The structure-losing behavior is even worse for cases like diff((x+1/x)^5,x); factor doesn't recover it.

[rtoy]

Imported from SourceForge on 2024-07-09 11:15:20 Created by rtoy on 2024-01-11 16:35:28 Original: https://sourceforge.net/p/maxima/bugs/4207/#0dbe/937d

---

Oh. That probably won't work because we're only differentiating once, and the result is pretty simple already. It's the combination of all the derivatives that is messy but simplifies nicely.