Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

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

2 Answers 2

up vote 6 down vote accepted

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.

share|improve this answer
    
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) 
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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