Re: How to compute functions defined by recurrences?

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Mon, Sep 23, 2024 at 05:05:54PM +0200, kevin lucas wrote:
> In some recent work I found myself having to investigate the asymptotics of
> the Debye polynomials, which I will here take as *defined* by the following
> integrodifferential recurrence:
> 
> $D_{n+1}(x)= \frac{1}{2}x^2(1-x)^2D'_n(x) +
> \frac{1}{8}\int_0^x(1-5t^2)D_n(t)dt$
> with $D_0(x)=1$

You can compute them with 

D(n)=
{
  my(d=1);
  for(i=1,n,d=x^2*(1-x)^2*(d')/2 + intformal((1-5*x^2)*d)/8);
  d
}

If you want the low-degree monomial only, you can do

Ds(n,l)=
{
  my(d=1+O(x^(l+1));
  for(i=1,n,d=x^2*(1-x)^2*(d')/2 + intformal((1-5*x^2)*d)/8);
  d
}

If you want the high degree term it is possible but requires
some manual computation (PARI does not know how to integrate 
x^n when n is a parameter).

Consider the operator from Q(x) to Q(x): 
F: $f \mapsto \frac{1}{2}x^2(1-x)^2f' + \frac{1}{8}\int_0^x(1-5t^2)fdt$

This linear operator has the property that 
if deg f > deg g, then deg F(f) > deg F(g),
which allows you to compute a formula for the leading term just by computing
the leading monomial of F(x^n) which is (n/2-5/(n+3)/8)*x^(3+n).

So we find that deg D_n = 3*n and that the

lc(D(n+1))*x^(3*(n+1))= lc(F(D(n)) = (3*n/2-5/(3*n+3)/8 * lc(D(n))

so we find that 

lc(D(n)) = prod(k=0, n-1 ,3*k/2-5/(3*k+3)/8)

So for example lc(D(1000)) is

? prod(n=0,999,3*n/2-5/(3*n+3)/8,1.)
%115 = -5.9847352224227929579864664579708515576E2739

You could use the same method for the coefficients of degree 3*n-1, 3*n-2 etc.

Cheers,
Bill