Double Sequences solving method

I am searching some information on what I think is called a real double sequence. What I call double sequence is a function from $\mathbb{N} \times \mathbb{N} \to \mathbb{R}$. By opposition to this, I call simple sequence (just a regular sequence), a function $\mathbb{N} \to \mathbb{R}$.

Recursion relation

I learnt that for a simple sequence, there are methods for solving recursion definitions of a sequence. For example, the simple function $u(n)$ defined as $\forall n\in\mathbb{N}, u(n+1)= u(n)+ u(n-1)$, we solve the characteristic polynomial associated with this recursion relation, and we find that the space of solution $S$ is the span of the Fibonacci numbers $F_n$ and Lucas numbers $L_n$ : $S\in Span(F_n,L_n)=Vect(F_n,L_n)$.

My problem

However, I have some interests in double sequences $u(i,j)$ defined by a recursion relation so that $\forall (i,j)\in \mathbb{N}^2,u(i,j) = f_1(i,j) \cdot u(i,j-1) + f_2(i,j) \cdot u(i-1,j-1)$, with the first term $u(0,0)$ being equal to a real constant, let's say $u(0,0) = \alpha$; and that $\forall (i,j)\in \mathbb{Z}^2, (i<0 \vee j<0),u(i,j) = 0$ so the terms with negative index in the recursion relation cancels out nicely.

A change of variable?

Someone suggested me to perform a change of variable so that we get a double sequence $v(i,j)$ with the form $\forall (i,j)\in \mathbb{N}^2,v(i,j) = g_1(i,j) \cdot v(i,j-1) + g_2(i,j) \cdot v(i-1,j)$. So, it makes sense for me because it seems so much easier to compute.

My questions

My questions are:

• Is there a general method for solving $ \forall (i,j)\in \mathbb{N}^2.$ $$ v(i,j) = g_1(i,j) \cdot v(i,j-1) + g_2(i,j) \cdot v(i-1,j) ?$$
• Is there a method if $g_1$ and $g_2$ are polynomials of second order of $i$ and $j$ : $g_1(i,j) = a + bi + ci^2 + dj + dj^2 + eij + fi^2j + gi^2j^2 + \dots $ ?
• Is there a matrix representation of this system that could be used to compute this?

Edit due to Yuval Filmus' comment

Edit: The comment turned out to be a little bit helpful with the user named Yuval Filmus telling me about generating functions and linking this book. Basically the idea of this clue is to replace the study of $u(i,j)$ by the study of $\sum \sum u(i,j) x^i y^j$; then, you replace the terms with partial derivatives as described on the book, and you turn your recursion problem in a PDE problem. However, I the problem is that the functions $g_1$ and $g_2$ are polynomials of order 2, and I don't know how to make it appear. So I would like to know if someone could explicitly turn my problem in a PDE problem then solve it? Then, once the function $\sum \sum u(i,j) x^i y^j$ has been found, you just do the do product with the family of functions $x^i$ and $y^j$ to get $u(i,j)$.

Side note: I didn't know about this method to solve recursion problem, but I knew the other way around: replace a PDE with $\sum \sum u(i,j) x^i y^j$ then identify the recursion relation, then solve it (assuming it was an easy one); then try to recognize it with the well-known Taylor expansions.

Here are two papers which might be relevant. The first one provides an analysis of linear recurrence relations in the multi-variate case based upon generating functions. The following abstract indicates that the geometry in more than one dimension makes things much more difficult than in the one-dimensional case, but also very interesting.

• Linear recurrences with constant coefficients: the multivariate case by Mireille Bousquet-Mélou and Marko Petkovšek

  Abstract: While in the univariate case solutions of linear recurrences with constant coeffcients have rational generating functions, we show that the multivariate case is much richer: even though initial conditions have rational generating functions, the corresponding solutions can have generating functions which are algebraic but not rational, D-finite but not algebraic, and even non-D-finite.

In section 4 the different types of solutions are analysed and criterias are stated from which the type of solution can be derived. Examples in the bivariate case are provided like (generalised) Dyck paths, Knight walks and some other types of paths.

The second recommendation provides an interesting analysis of a subclass of bivariate linear recurrences, which was stated as research problem in Concrete Mathematics by R.L. Graham, D.E. Knuth and O. Patashnik.

• Bivariate Generating Functions for a Class of Linear Recurrences: General Structure by J. Fernando Barbero G., Jesús Salas and Eduardo J. S. Villaseñor.

An interesting aspect is that this class of linear recurrences covers many well-known numbers, like Lah numbers, Stirling numbers of both kinds, Eulerian numbers and many more.

In fact there is a total of $23$ different integer sequences stated which are classified according to four different types of partial differential equations (PDEs). Each of these PDEs is analysed providing this way instructive examples how to cope with generating functions, resp. with their underlying recurrence relations.

Application: Knight's walk

Here we go through parts of the first paper, look at the seemingly innocent bivariate recurrence relation of the knight's walk and work out how to classify such relations.

We start with the recurrence relation stated in the multivariate case and the bivariate specialisation for the knight's walk \begin{align*} a_{\mathbf{n}}&= \begin{cases} \sum_{\mathbf{h}\in H}c_\mathbf{h}a_{\mathbf{n}+\mathbf{h}}&\qquad\qquad\quad \mathbf{n}\geq \mathbf{s}\\ \phi(\mathbf{n})&\qquad\qquad\quad \mathbf{n}\geq \mathbf{0},\mathbf{n\not \geq s}\tag{1} \end{cases}\\ \color{blue}{a_{m,n}}&= \begin{cases} \color{blue}{a_{m+1,n-2}+a_{m-2,n+1}}&\qquad \color{blue}{m,n\geq 2}\\ \color{blue}{1}&\qquad \color{blue}{m,n\geq 0, (m,n)\not\geq (2,2)}\tag{2} \end{cases} \end{align*}

The bivariate recurrence relation (2) of knight's walk is a specialisation of the multivariate case (1) with
\begin{align*} H=\{(1,-2),(-2,1)\}\qquad \text{and}\qquad s=(2,2) \end{align*} The inequality relation $(m,n)\not\geq (2,2)$ is a component-wise order relation meaning that not both coordinates are allowed to be greater or equal $2$.

We define \begin{align*} F(\mathbf{x})&=\sum_{\mathbf{n}\geq \mathbf{0}}a_{\mathbf{n}}\mathbf{x}^\mathbf{n}\qquad\qquad \color{blue}{F(x,y)=\sum_{{m\geq 0}\atop{n\geq 0}}a_{m,n}x^my^n} \end{align*} and consider $F_s(\mathbf{x})=\sum_{\mathbf{n}\geq \mathbf{s}}x^{\mathbf{n}-\mathbf{s}}$.

We obtain

\begin{align*} \color{blue}{F_s(x,y)}&\color{blue}{=\sum_{{m\geq 2}\atop{n\geq 2}}a_{m,n}x^{m-2}y^{n-2}}\\ &=\sum_{{m\geq 2}\atop{n\geq 2}}\left(a_{m+1,n-2}+a_{m-2,n+1}\right)x^{m-2}y^{n-2}\tag{3}\\ &=\frac{y^2}{x}\sum_{{m\geq 3}\atop{n\geq 0}}a_{m,n}x^{m-2}y^{n-2}+\frac{x^2}{y}\sum_{{m\geq 0}\atop{n\geq 3}}a_{m,n}x^{m-2}y^{n-2}\tag{4}\\ &=\frac{y^2}{x}\left[F_s(x,y)+\sum_{m\geq 3}\left(a_{m,0}x^{m-2}y^{-2}+a_{m,1}x^{m-2}y^{-1}\right)\right.\\ &\qquad\qquad\quad\left.-\sum_{n\geq 2}a_{2,n}x^0y^{n-2}\right]\tag{5}\\ &\qquad+\frac{x^2}{y}\left[F_s(x,y)+\sum_{n\geq 3}\left(a_{0,n}x^{-2}y^{n-2}+a_{1,n}x^{-1}y^{n-2}\right)\right.\\ &\qquad\qquad\quad\left.-\sum_{m\geq 2}a_{m,2}x^{m-2}y^0\right]\\ &=\frac{y^2}{x}F_s(x,y)+(1+y)\sum_{m\geq 3}x^{m-3}-\frac{1}{x}\sum_{n\geq 2}a_{2,n}y^{n-2}\tag{6}\\ &\qquad+\frac{x^2}{y}F_s(x,y)+(1+x)\sum_{n\geq 3}y^{n-3}-\frac{1}{y}\sum_{n\geq 2}a_{m,2}x^{m-2}\\ &\,\,\color{blue}{=\left(\frac{x^2}{y}+\frac{y^2}{x}\right)F_s(x,y)+\left(\frac{1+y}{1-x}+\frac{1+x}{1-y}\right)}\\ &\qquad\color{blue}{-\left(\frac{1}{x}\sum_{n\geq 2}a_{2,n}y^n+\frac{1}{y}\sum_{m\geq 2}a_{m,2}x^{m}\right)}\tag{7}\\ \end{align*}

Comment:

• In (3) we use the recurrence relation (2).

• In (4) we shift the indices to get coefficients $a_{m,n}$.

• In (5) we separate $F_s(x,y)$.

• In (6) we use that boundary values at lines $x=0,1$ and $y=0,1$ are $1$.

• In (7) we do a geometric series expansion and collect terms.

We take the points $(1,-2), (-2,1)$ from $H$ and derive a point $$\mathbf{p}=(p_1,p_2), p_i:=\max\{h_i,1\}, i=1,2$$ called apex, which is

\begin{align*} \color{blue}{P=(1,1)}\tag{8} \end{align*}

We multiply (7) with $x^{\mathbf{p}}=xy$ and obtain

\begin{align*} \color{blue}{F_s(x,y)}&\color{blue}{=\frac{K(x,y)-U(x,y)}{Q(x,y)}}\tag{9}\\ Q(x,y)&=xy-x^3-y^3\\ K(x,y)&=xy\left(\frac{1+y}{1-x}+\frac{1+x}{1-y}\right)\\ U(x,y)&=\sum_{m\geq 2}a_{m,2}x^{m+1}+\sum_{n\geq 2}a_{2,n}y^{n+1} \end{align*}

The function $Q(x,y)$ in (9) is called the kernel function and is completely known as well as $K(x,y)$, the known function. The knowledge of $F_s(x,y)$ depends solely on the unknown function $U(x,y)$.

If $F_s(x,y)$ is known, we also know $F(x,y)$, since

\begin{align*} F(x,y)=x^2y^2\left(F_s(x,y)+\sum_{{m,n\geq 0}\atop{(m,n)\not\geq (2,2)}}a_{m,n}x^my^n\right) \end{align*}

We are now well prepared for a classification of this recurrence relation and linear recurrence relations in general. Interestingly, the apex $P=(1,1)$ from (8) plays a key role.

The following is valid

• (from Theorem 12): Assume the apex $P=(p_1,p_2)$ of $H$ is $(0,0)$. Then the generating function $F_s(x,y)$ of the unique solution of (2) is rational if and only if the known initial function $K(x,y)$ itself is rational.

• (from Theorem 13): Assume the apex $P=(p_1,p_2)$ of $H$ has at most one positive coordinate. Then the generating function $F_s(x,y)$ of the unique solution (2) is algebraic if and only if the known initial function $K(x)$ itself is algebraic.

Since the apex $P=(1,1)$ has more than one positive coordinates, we conclude from these theorems that the knight's walk has neither a rational nor an algebraic generating function. In fact the knight's walk has not even a $D$-finite generating function. This is shown at the end of the paper together with a complete solution of the recurrence relation.

The key word in this question is solve. What does it mean to "solve" a recursion relation? Looking at a few simple examples, it is clear that this is like asking to find integrals of general functions "in closed form". The answer is that you can not do this except for some special cases.

Considering just a simple one dimensional recursion relation $\ u(n+1) = f(u(n)),\ $ with initial value $\ u(0)=z,\ $already you can see that although you can compute any value $\ u(n)\ $ by starting from $\ u(0),\ $ in general, there is no better way to find the value. Thus, for your first question about general methods for solving recursion relations even in one dimension or higher, the answer is no. For your second question, even restricting the functions $\ g_1()\ $ and $\ g_2()\ $ to be polynomials at most of second degree in $\ i,j\ $ the answer is still no.

For your "double sequences" a few preliminary remarks. Consider the simplest example of such sequences which is the binomial coefficients, also known as Pascal's triangle. In one version they are defined by $\ u(n,k) := (n+k)!/(n!k!)\ $ and are the entries of a symmetric infinite square matrix with recursion relation $\ u(n,k) = u(n-1,k) + u(n,k-1).\ $ In a second version they are defined by $\ u(n,k) := n!/(k!(n-k)!)\ $ and are the entries of a symmetric infinite triangular array with recursion relation $\ u(n,k) = u(n-1,k) + u(n-1,k-1).\ $ In either version, the indexing could either begin with $\ n=k=0\ $ or with $\ n=k=1.\ $ These are merely matters of personal choice and convenience.

Let us try solving a few simple examples of the recursion $$ v(i,j) = g_1(i,j) \cdot v(i,j-1) + g_2(i,j) \cdot v(i-1,j) $$ with initial value $\ v(1,1) = 1.\ $

First, let $\ g_1(i,j) = 1, \ g_2(i,j) = i.\ $ The array of numbers is OEIS sequence A132159. Explicitly, the solution is $$ v(i,j) = i\ (i+j-2)!/(j-1)!. $$ Using the information in the OEIS entry, the exponential generating function (egf) is given by $$ \frac{\exp(y z)}{(1 - xz)^2} = \sum_{n=0}^\infty \frac{z^n}{n!}\: \Big(\sum_{i+j=n+2} v(i,j)\ x^i y^j\Big). $$

Next, let $\ g_1(i,j) = j, \ g_2(i,j) = i.\ $ The array of numbers is OEIS sequence A091441. Explicitly, the solution is $$ v(i,j) = i\ j\ (i+j-2)!. $$ Using the information in the OEIS entry, the egf is given by $$ \frac1{((1 - xz)(1-yz))^2} = \sum_{n=0}^\infty \frac{z^n}{n!}\: \Big(\sum_{i+j=n+2} v(i,j)\ x^i y^j\Big). $$

Next, let $\ g_1(i,j) = i, \ g_2(i,j) = 1.\ $ The array of numbers is OEIS sequence A106800. Explicitly, the solution is $$ v(i,j) = S2(i+j-1,i) $$ where $S2()$ is the Stirling numbers of the second kind. Using the information in the OEIS entry, the egf is given by $$ \exp((\exp(y\ z)-1)\ x/y)-1 = \sum_{n=1}^\infty \frac{z^n}{n!}\: \Big(\sum_{i+j=n} v(i,j)\ x^i y^j\Big). $$

Next, let $\ g_1(i,j) = i, \ g_2(i,j) = j.\ $ The array of numbers is OEIS sequence A008292. These are the Eulerian numbers. Using the information in the OEIS entry, the egf is given by $$ \frac{x-y}{1-y/x\ \exp(xz-yz)}-x = \sum_{n=1}^\infty \frac{z^n}{n!}\: \Big(\sum_{i+j=n+1} v(i,j)\ x^i y^j\Big). $$

There are several other simple examples of this kind that I have worked on so far, but my usual method is to compute the first few numbers and look them up in the OEIS and use the information there to find an explicit expression for the numbers and/or an exponential generating function for them. This only works in certain special cases, of course.

For example, let $\ g_1(i,j) = 1, \ g_2(i,j) = i^2.\ $ The array of numbers is currently not in the OEIS, but despite this, the explicit expression is $\ v(i,j) = i\ i!\ (i+j-2)!/(j-1)!\ $ but so far I have not found an expression for the egf.

My experience strongly suggests that quadratic functions $\ g_1(),\ g_2()\ $ are much harder to solve recursions for. I don't know of any example of PDE that are associated with these quadratic functions, and even then, the PDE may not be any easier to solve for.