Really hard two-dimensial recurrence relations of form p[n,m]=p[n,m-1]+p[n-1,m]+p[n-1,m-1]*(n-1)

Question

I have a two-dimensional recurrence equation, help me solve this:

p[n,m]=p[n,m-1]+p[n-1,m]+p[n-1,m-1]*(n-1)

p[n,0]=1 

p[0,m]=0

p[0,0]=0

I generated these numbers for 1<=n,m<=6:

n row, m column

1 1 1 1 1 1

3 5 7 9 11 13

6 17 34 57 86 121

10 45 130 289 546 925

15 100 410 1219 2921 6030

21 196 1106 4375 13391 34026

Firstly I saw, that p[n,1] = n*(n+1)/2

Next, fix n = 2, look for the differences between p[n,i] and p[n,i-1].

They are all equals 2 = 2! (remember that)

Now, fix n = 3, also look for the differences between p[n,i] and p[n,i-1]

We have 11, 16, 23, 29. Okay so now look for the differences between differences :)

They are all equals 6 = 3!

Now, fix n = 4, also (hah) look for the differences between p[n,i] and p[n,i-1]

We have 35, 85, 159, 257. Look for the differences between differences.

We have 50, 74, 98. Also look for the differences between differences.

They are all equals 24 = 4!

Now, fix n = 5, also (hah) look for the differences between p[n,i] and p[n,i-1]

85, 310, 809, 1702 ->

225, 499, 893 ->

274, 394 ->

120 = 5!

And so on...

That's all for now :(

updated: I found oeis sequence which is very similar to mine!

Answer 1

I suspect the solution to this difference equation is very strongly (factorially?) divergent, so calculating values for n,m as large as 10^5 is going to be challenging.

One can obviously (and inefficiently!) compute p(n,m) by simple recurrence as follows (in Python):

import numpy
n_max, m_max = 10, 10
p = numpy.zeros((n_max+1, m_max+1), dtype='int64')
p[:,0] = 1
p[0,:] = 0
p[0,0] = 0
for n in range(1, n_max+1):
    for m in range(1, m_max+1):
        p[n,m] = p[n, m-1] + p[n-1, m] + p[n-1, m-1] * (n-1)

This gives the following result for p(n,m):

[[         0          0          0          0          0          0
           0          0          0          0          0]
 [         1          1          1          1          1          1
           1          1          1          1          1]
 [         1          3          5          7          9         11
          13         15         17         19         21]
 [         1          6         17         34         57         86
         121        162        209        262        321]
 [         1         10         45        130        289        546
         925       1450       2145       3034       4141]
 [         1         15        100        410       1219       2921
        6030      11180      19125      30739      47016]
 [         1         21        196       1106       4375      13391
       34026      75356     150381     276745     477456]
 [         1         28        350       2632      13643      53284
      167656     447168    1049685    2228716    4366642]
 [         1         36        582       5664      37731     186516
      727160    2347920    6527781   16104292   36071946]
 [         1         45        915      11235      94278     582642
     2801930   10967130   36278271  104604811  269511093]
 [         1         55       1375      20845     216238    1647382
     9693090   45877590  180860031  611969281 1822923673]]

which already contains some sizeable values even for n=m=10. Extending that computation to n=m=100, and using floating-point arithmetic, indicates that p(100,100) may be as large as 5x10^172.

Using a generating function

,

I believe you can convert your 2D difference equation into something like

,

which perhaps might help your analysis. However, as an illustrative comparison, one could consider a difference equation of the form

which can be converted to the following differential equation for the generating function:

which has a solution of the form:

Clearly, such a generating function would have a very poorly behaved Taylor expansion near .