# Implementation of a recurrence relation for the polynomials appearing in the large order asymptotics of the Bessel functions

I would like to implement the recurrence relation for the polynomials $U_n(x)$ appearing in the large order asymptotics of the Bessel functions. The recurrence in question is $$U_{n+1}(x)=\frac{1}{2}x^2(1-x^2)U_n^{\prime}(x)+\frac{1}{8}\int_{0} ^x(1-5t^2)U_n(t)dt,$$ with $U_0(x)=1$. I saw implementations for recursions for polynomials, but this one involves integration and differentiation. Any idea? Thank you very much in advance!

-
Related problem concerning Hermite polynomials: Why does Expand not work within a function? –  Artes Jul 25 '13 at 21:00

I tried to find a faster way to generate these coefficients, compared to memoization. What I came up with is this:

With[ { n = 3},
Last[
coefficientList =
NestList[Simplify[1/2 x^2 (1 - x^2) D[#, x] +
1/8 Subtract @@ ({#, # /. x -> 0} &@
Integrate[(1 - 5 x^2) #, x])] &, 1, n]]
]

(*
==> 1/
8 ((3 x^3)/128 - (289 x^5)/1920 + (385 x^7)/1152 - (1925 x^9)/
10368) +
1/2 x^2 (1 - x^2) (-(1/8) x^3 (1 - 5 x^2) - 5/8 x^3 (1 - x^2) +
1/8 x (1 - 5 x^2) (1 - x^2) + 1/8 (x/8 - (5 x^3)/6 + (25 x^5)/24))
*)


Here, the Last at the top is only there to suppress the output of the entire list of coefficients. The number of coefficients is n = 3. There are two things that make this significantly faster than the recursive definition by memoization:

First, NestList is a built-in function for doing recursions. Second, the Integrate is faster when done as an indefinite integral instead of a definite integral. So I do the former and substitute the integration boundaries afterwards.

This is about an order of magnitude faster than memoization for n = 4.

Edit

To get even better performance at large n, it's also good to add Simplify at each step:

With[{n = 10},
Last[coefficientList =
NestList[
Simplify[
1/2 x^2 (1 - x^2) D[#, x] +
1/8 Subtract @@ ({#, # /. x -> 0} &@
Integrate[(1 - 5 x^2) #, x])] &, 1, n]]]

(*
==> (x^10 (9745329584487361980740625 -
1230031256571145165088463750 x^2 +
27299183373230345667273718125 x^4 -
246750339886026017414509498824 x^6 +
1177120360439828012193658602930 x^8 -
3327704366990695147540934069220 x^10 +
5876803711285273203043452095250 x^12 -
6564241639632418015173104205000 x^14 +
4513386761946134740461797128125 x^16 -
1745632061522350031610173343750 x^18 +
290938676920391671935028890625 x^20))/88580102706155225088000
*)


This makes the integrals easier to calculate at larger n because the number of terms is reduced.

-
I guess Expand would be faster and quite sufficient (+1). –  Artes Jul 25 '13 at 21:17
@Artes Indeed, Simplify is a little bit slower then Expand here. BTW, I had also upvoted your hermite answer using Expand, as I just saw... –  Jens Jul 25 '13 at 21:24
Thank you very much! One more question: how can I evaluate the resulting polynomial at a given complex number? –  Gary Jul 26 '13 at 6:19
Ok, I've found the answer myself. –  Gary Jul 26 '13 at 6:23

You can use memoization :

u[0, x_] = 1;
u[n_, x_] := u[n, x] = 1/2 x^2 (1 - x^2) D[u[n - 1, x], x] +
1/8 Integrate[(1 - 5 t^2) u[n - 1, t], {t, 0, x}]


Check :

u[2, x]

1/16 x^2 (1 - 5 x^2) (1 - x^2) + 1/8 (x^2/16 - (5 x^4)/24 + (25 x^6)/144)

-