RSolve two variables with multiple boundary conditions

Question

In short, I want to solve a recurrence equation in two variables with multiple boundary conditions.

In more detail:

For instance, suppose I want to compute

$a[0,j]= 1 \quad \forall j$

$a[i,0]= 1 \quad \forall i$

$a[i,j] = a[i-1,j]+a[i-1,j-1]+a[i,j-1]\quad \forall i,j \geq 1$

My approach:

Clear["Global`*"]
RSolve[{a[i, j] == a[i - 1, j] + a[i - 1, j - 1] + a[i, j - 1], a[0, j] == 1, a[i, 0] == 1}, a[i, j], {i, j}]

Just yields the same expression as output -- without any error message or result:

RSolve[{a[i, j] == a[-1 + i, -1 + j] + a[-1 + i, j] + a[i, -1 + j], a[0, j] == 1, a[i, 0] == 1}, a[i, j], {i, j}]

I presume the error is because I somehow need to encode the case distinction (for $i,j \geq 1$, apply the recursive equation, else the constants), but I could not find any example code in the documentation on how to do this. In general, in the documentation, there is never a case distinction, so mathematica seems to ''find out'' which equation to use.

My attempt to fix this was a piecewise function:

RSolve[a[i, j] == Piecewise[{{a[i - 1, j], i > 0}, {1, i == 0} } ], a[i, j], {i, j}]

But also did not work.

Debugging the code further, I have the weird behaviour that this very basic version of the recurrence equation works:

RSolve[{a[i, j] == a[i - 1, j], a[0, j] == 1}, a[i, j], {i, j}]

{{a[i, j] -> 1}}

but as soon as I add the second base case:

RSolve[{a[i, j] == a[i - 1, j], a[0, j] == 1, a[i, 0] == 1},  a[i, j], {i, j}]

I get again the strange output

RSolve[{a[i, j] == a[-1 + i, j], a[0, j] == 1, a[i, 0] == 1}, a[i, j], {i, j}]

without any error message.

My questions would be:

(1) Does anyone know how to solve this problem?

(2) Going further, if we have another case distinction, e.g. the recurrence would change for $i > 10 $ to $ a[i,j] = 100*a[i-1,j]+ a[i-1,j-1] + a[i,j-1]$, is there a way to encode that case distinction?

Answer 1

Not a solution with RSolve, but get recurrenceTable with Reduce and find " Delannoy numbers d(i,j) " as generating formula.

eqs = {a[i, j] == a[i - 1, j] + a[i - 1, j - 1] + a[i, j - 1], 
       a[0, j] == 1, a[i, 0] == 1}

vars = Flatten[mat = Table[a[i, j], {i, 0, 6}, {j, 0, 6}]];

tab = Flatten@Table[eqs, {i, 0, 6}, {j, 0, 6}] // Union

red = Reduce[tab, vars, Reals]

sol = Solve[red]

(mat1 = mat /. sol[[1]]) // MatrixForm

(*   MatrixForm[{{1, 1, 1, 1, 1, 1, 1}, 
                 {1, 3, 5, 7, 9, 11, 13}, 
                 {1, 5, 13, 25, 41, 61, 85}, 
                 {1, 7, 25, 63, 129, 231, 377}, 
                 {1, 9, 41, 129, 321, 681, 1289}, 
                 {1, 11, 61, 231, 681, 1683, 3653}, 
                 {1, 13, 85, 377, 1289, 3653, 8989}}]   *)

And @@ (tab /. sol[[1]])

(*   True   *)

With the help of https://oeis.org/search?q=1%2C13%2C61%2C129&sort=&language=english&go=Search you find, the results are the " Delannoy numbers d(i,j) ".

d[n_, k_] := Binomial[n + k, k]*Hypergeometric2F1[-k, -n, -n - k, -1];

Table[d[n, k], {n, 0, 6}, {k, 0, n}] // MatrixForm

Since solution matrix is symmetric, get it with

(mat2 = Table[a[i, j], {i, 0, 6}, {j, 0, 6}] /. 
    a[i_, j_] -> d[i, j]) // MatrixForm

mat1 == mat2
(*   True   *)