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?

  • 4
    $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Jul 22 at 9:34
  • $\begingroup$ Thank you very much! That seems to be working and giving the sequences' terms. $\endgroup$
    – 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


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<n\leq 5$ 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}.

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.