# Implementation of a complex recurrence relation for polynomials

I would like to implement the recurrence relation for the polynomials $U_n(x)$ that appear 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 of recursions for polynomials, but this one involves integration and differentiation. Any ideas?

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). Jul 25, 2013 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, 2013 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, 2013 at 6:19
• Ok, I've found the answer myself.
– Gary
Jul 26, 2013 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)


One more possibility would be to directly manipulate lists of coefficients of polynomials. To do that, we need to define a few required operations:

SetAttributes[{add, mult}, Orderless];
mult[{0}, c2_?VectorQ] := {0};
mult[c1_?VectorQ, c2_?VectorQ] := ListConvolve[c1, c2, {1, -1}, 0];

diff[{_}] := {0};
diff[c_?VectorQ] := Rest[c] Range[Length[c] - 1];
int[c_?VectorQ] := Prepend[c/Range[Length[c]], 0];


and with those in hand,

NestList[add[mult[{0, 0, 1/2, 0, -1/2}, diff[#]], int[mult[{1/8, 0, -5/8}, #]]] &, {1}, 5]
{{1}, {0, 1/8, 0, -5/24}, {0, 0, 9/128, 0, -77/192, 0, 385/1152},
{0, 0, 0, 75/1024, 0, -4563/5120, 0, 17017/9216, 0, -85085/82944},
{0, 0, 0, 0, 3675/32768, 0, -96833/40960, 0, 144001/16384, 0, -7436429/663552,
0, 37182145/7962624},
{0, 0, 0, 0, 0, 59535/262144, 0, -67608983/9175040, 0, 250881631/5898240, 0,
-108313205/1179648, 0, 5391411025/63700992, 0, -5391411025/191102976}}


Use FromDigits[Reverse[coeff], x] to reconstitute the polynomials.