# Implementing solution of a system of recurrence relations

Let us consider a system of recurrence relations such as

$$a_{n-1} = ( \lambda_1 + n \lambda_2 ) a_n + \lambda_3 b_n$$

$$b_{n-1} = ( \lambda_4 + n \lambda_2 ) b_n + \lambda_5 a_n$$

subject to the initial conditions say, $$a_0 = k_1$$ and $$b_0 = k_2$$. So, what I have done so far is that I reexpressed the system as

$$v_n = A^T v_{n-1}$$

or equivalently,

$$v_n = (A^T)^n v_0$$,

where $$v_n = (a_n, b_n)^T$$, and $$A$$ is the coefficient matrix including $$n$$.

How do I implement such a system in Mathematica to get the sequences $$a_n$$ and $$b_n$$?

Also, is there a function, in Mathematica, similar to RSolve for solving such systems?

• Unfortunately RSolve[Join[Thread[{a[n - 1], b[n - 1]} == {{λ1 + n λ2, λ3}, {λ5, λ4 + n λ2}} . {a[n], b[n]}], Thread[{a[0], b[0]} == {k1, k2}]], {a, b}, n] doesn't work. However, you can always do RecurrenceTable[Join[Thread[{a[n - 1], b[n - 1]} == {{λ1 + n λ2, λ3}, {λ5, λ4 + n λ2}} . {a[n], b[n]}], Thread[{a[0], b[0]} == {k1, k2}]], {a, b}, {n, 0, 10}], for example. – J. M.'s torpor Jul 22 at 9:34
• Thank you very much! That seems to be working and giving the sequences' terms. – lynx Jul 22 at 11:12

Because the coefficient matrix commutes with itself for all $$n$$ - you can solve this as a scalar equation then generalize the result. This is not exactly general and a little inelegant, but I can't seem to get Mathematica to do matrix recurrance relations natively. It would be nice if this were supported for general operators.

Anyway we start by noting that

$$a_{n-1} = ( \lambda_1 + n \lambda_2 ) a_n + \lambda_3 b_n \quad \text{and} \quad b_{n-1} = ( \lambda_4 + n \lambda_2 ) b_n + \lambda_5 a_n$$

can be written as $$A_{n-1}=\tilde B_n A_n=(B+n\lambda_2 I)A_n\quad \text{with}\quad B= \begin{pmatrix} \lambda _1 & \lambda _3 \\ \lambda _5 & \lambda _4 \end{pmatrix}$$

Pretending this is a scalar equation lets us solve

scalarSol=A[n]/.Simplify[First@RSolve[{A[n-1]==(B+n λ2)A[n],A[0]==k},A[n],n],{n∈Integers}]


which returns $$\dfrac{k \lambda _2^{1-n}}{\left(B+\lambda _2\right) \left(\frac{B}{\lambda _2}+2\right)_{n-1}}$$ where $$\left(\cdot\right)_n$$ is the Pochhammer symbol. Mathematica struggles to simplify this but can confirm it's equivalent to $$\dfrac{k}{\lambda _2^{n}\left(\frac{B}{\lambda _2}+1\right)_n}$$, i.e.

FullSimplify[scalarSol == k λ2^-n/Pochhammer[1 + B/λ2, n]] gives True.

From this, we turn scalarSol into a matrix expression. $$B$$ is the matrix given above, $$k=A_0=\begin{pmatrix} a_0 \\b_0\end{pmatrix}=\begin{pmatrix} k_1 \\k_2\end{pmatrix}$$ is the initial condition and $$1$$ is the identity matrix. In general this transformation is not well defined but, again, since $$(B+n\lambda_2 I)$$ commutes with $$(B+m\lambda_2 I)$$ for all $$n$$ and $$m$$ we don't have to worry about most of the subtleties that typically arise. The only question would be if $$k$$ should be multiplied on the left or right, and it should obviously be the rightmost quantity here. Therefore we have

sol[n_]=FullSimplify[With[{B={{λ1,λ3},{λ5,λ4}}},
Inverse[MatrixFunction[Pochhammer[#,n]&,B/λ2+IdentityMatrix[2]]].{k1,k2}]/λ2^n]


which is an explicit, analytic solution $$\begin{pmatrix} a_n\\b_n \end{pmatrix} = \begin{pmatrix} \frac{\left(\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5}+\lambda _1-\lambda _4\right) k_1+2 k_2 \lambda _3}{\left((\lambda _1+2 \lambda _2+\lambda _4+\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5})/(2 \lambda _2)\right){}_n} + \frac{k_1 \left(\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5}-\lambda _1+\lambda _4\right)-2 k_2 \lambda _3}{\left((\lambda _1+2 \lambda _2+\lambda _4-\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5})/(2 \lambda _2)\right){}_n} \\ \frac{k_2 \left(\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5}+\lambda _1-\lambda _4\right)-2 k_1 \lambda _5}{\left((\lambda _1+2 \lambda _2+\lambda _4-\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5})/(2 \lambda _2)\right){}_n} + \frac{\left(\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5}-\lambda _1+\lambda _4\right) k_2+2 k_1 \lambda _5}{\left((\lambda _1+2 \lambda _2+\lambda _4+\sqrt{\left(\lambda _1-\lambda _4\right){}^2+4 \lambda _3 \lambda _5})/(2 \lambda _2)\right){}_n} \end{pmatrix}$$

We can check this matches for $$0 and indeed

check[n1_]:=First@RecurrenceTable[Join[Thread[{a[n-1],b[n-1]}=={{λ1+n λ2,λ3},{λ5,λ4+n λ2}}.{a[n],b[n]}],Thread[{a[0],b[0]}=={k1,k2}]],{a,b},{n,n1,n1}]
Table[check[m] == sol[m] // FullSimplify, {m, 1, 5}]


returns {True,True,True,True,True}.

• Thank you very much! This is quite useful indeed. – lynx Jul 23 at 7:52