Re: How to compute functions defined by recurrences?

k is -0.120526166805581503451769239491718240700000837213342230433749796337295780027993182203117298864948031052323469664... thanks to PARI's fast prodnumrat.

I started by computing the leading coefficients and saw they were growing quickly, so I thought to look at their first quotients. These seemed to be growing linearly, so I looked at their differences and it seemed to settle down pretty quickly to 1.5 so I was fairly sure the main term would be something like 1.5k + a for some a. Solving for a with as large a term as I could quickly compute I saw it was suspiciously close to -1.5, so the total looked like 1.5 * 3 * 4.5 * ... * 1.5(n-1) which is just (n-1)! times 1.5^(n-1). I chose to write 1.5^n because it was neater and absorb the extra 1.5 into the constant.

On Mon, Sep 23, 2024 at 2:23 PM kevin lucas <lucaskevin296@gmail.com> wrote:

Thanks so much for your replies! I’m amazed how fast you got somewhere and still trying to digest all of it, but I think I have more than enough to continue with for now. To Charles, how did you guess the (n-1)! for the leading coefficient?

Also if there's a way to get enough digits of k to try to recognize it I’d be very interested.

On Mon 23. Sep 2024 at 19:09, Charles Greathouse <crgreathouse@gmail.com> wrote:

And indeed
-5/36 * prodnumrat(1-5/24/k/(3*(k-1)/2),2)

is about -0.1205, so my empirical k was (approximately) correct.

I recognize that this doesn't address the part of your question about computing the whole polynomial for large n. My code generates D_250(x) in about 8 seconds, which isn't bad considering the output is around 1 MB (a bit more printed, a bit less internal). I'm sure there are clever ways to do it better but if you need a bit output it's never going to be super fast.

On Mon, Sep 23, 2024 at 12:36 PM Charles Greathouse <crgreathouse@gmail.com> wrote:

I'm not afraid of inflicting bad code (or embarrassing myself) on the list, so I'll give

D(n)=if(n==0, return(1)); my(d=D(n-1)); x^2*(1-x)^2*(d')/2 + intformal((1-5*x^2)*d)/8

which I think is the naive transcription of the definition. Let's use f(n) as the leading coefficient of D_n(x). Then, speaking purely numerically, it looks like

f(n)/f(n-1) = 1.5n + C

where the constant looks suspiciously like -1.5. If that's right, you should have a factorial-like function as the leading coefficient, something like

f(n) = -k * (n-1)! * (3/2)^n

Numerically it seems to work pretty well with a constant k near 0.120.

OK, so how did we do? From the defining formula, it looks like

f(n) = prod(k=1,n, (3*k-3)/2-5/24/k)

which indeed has the (n-1)! * (3/2)^n but I missed the correction term.

On Mon, Sep 23, 2024 at 11:06 AM kevin lucas <lucaskevin296@gmail.com> 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$

These are polynomials of degree 3n. I mostly care about the leading term as n grows, but I'd also like to know how to compute these polynomials to high degree assuming only this recurrence. I have a horribly slow and ugly GP script which I am too embarrassed by to inflict on this list, and it doesn't go high enough to get the digits to guess at asymptotics.

More generally, I can't find from the usual sources how to accelerate computing families of *functions* given by recurrences. If for example all one had was this recurrence for Bernoulli polynomials, could one easily compute B_1000(x) (or at least an arbitrary coefficient) in PARI/GP?

As usual any help, comments or references are appreciated.