Bessel derivatives

Question

I tried to expand BesselJ[k,x] function into a Taylor series with Series command. Here both k and x are some functions of the expansion variable $\lambda$, so in the expansion, derivatives with respect to both k and x occur.

The problem is, whenever there is a term that is the derivative of both variables, Mathematica leaves it as (e.g.) Derivative[2, 1][BesselJ][0., 2.40483] and doesn't give a numerical value in the end.

First I thought it is because assigned values aren't exact -- since the above 2.40483 is the value of BesselJZero[0, 1]. However, it also doesn't give numerical or analytical result for the following easier calculation:

D[BesselJ[k, x], {k, 1}, {x, 1}] /. {k -> 1, x -> 1} // N
(* -> Derivative[1, 1][BesselJ][1., 1.] *)

But, when the order of the differentiation changes, it gives numerical results:

D[BesselJ[k, x], {x, 1}, {k, 1}] /. {k -> 0, x -> 1} // N
(* -> 1.22713 *)

It works with x -> BesselJZero[0, 1] as well.

First question: why is this the case?

Second (if it is possible): how can I handle it with Series command?

Answer 1

Short answer: Use FunctionExpand.

Longer answer: This indeed a bit odd at first glance. However, it seems to be more a property of the Bessel functions than of D. In any case, using FunctionExpand solves the problem:

D[BesselJ[k, x], k, x] /. {k -> 1, x -> 1} // N
FunctionExpand[%]

yields

0.160512

Now, that seems to work in a series as well:

Normal@Series[BesselJ[ Sin[x], Cos[x]], {x, 0, 3}] // N

gives

 0.7651976865579666 + 0.13863371520404605*x - 0.6542004860899648*x^2 + 
 0.16666666666666666*x^3*(2.5114247538954846 - 3.*Derivative[1, 1][BesselJ][0., 1.])

Applying FunctionExpand on it, resolves the remaining symbolic derivatives

FunctionExpand[%]
0.7651976865579666 + 0.13863371520404605*x - 0.6542004860899648*x^2 -  0.19499232275587164*x^3

Answer 2

You can get numerical series expansions by using the NumericalCalculus package, which offers the function NSeries. I'll illustrate this with some simple made-up functions on which the BesselJ arguments depend, as you describe. Using ϵ as the small parameter, I define the fifth-order series expansion as f[ϵ_] and compare it to the exact function.

To make the numerics work reliably, I have to specify a radius on which the function is sampled. I choose this small enough to exclude the potentially troublesome points where the Bessel arguments are zero.

<< NumericalCalculus`

f[ϵ_] = 
 Chop@Normal[
   NSeries[BesselJ[1 + ϵ^2, 
     Sin[Pi/4 + ϵ]], {ϵ, 0, 5}, Radius -> .5]]

(*
==> 0.331912 + 0.289531 ϵ - 0.681305 ϵ^2 - 
 0.0643912 ϵ^3 + 0.546373 ϵ^4 - 0.636954 ϵ^5
*)

Plot[{f[ϵ], 
  BesselJ[1 + ϵ^2, Sin[Pi/4 + ϵ]]}, {ϵ, 
  -.5, .5}, PlotStyle -> {Dashed, Red}]

This shows that the expansion worked as it should.