Implementation of a complex recurrence relation for polynomials

Question

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?

Answer 1

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.

Answer 2

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) 

Answer 3

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];
add[c1_?VectorQ, c2_?VectorQ] := Total[PadRight[{c1, c2}]];
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.