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?

  • 3
    $\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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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