6
$\begingroup$

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?

$\endgroup$
10
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I guess Expand would be faster and quite sufficient (+1). $\endgroup$ – Artes Jul 25 '13 at 21:17
  • $\begingroup$ @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... $\endgroup$ – Jens Jul 25 '13 at 21:24
  • $\begingroup$ Thank you very much! One more question: how can I evaluate the resulting polynomial at a given complex number? $\endgroup$ – Gary Jul 26 '13 at 6:19
  • $\begingroup$ Ok, I've found the answer myself. $\endgroup$ – Gary Jul 26 '13 at 6:23
7
$\begingroup$

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) 
$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.