How to calculate the equilibrium of a difference equation system of multiple functions?

Question

I'm trying to infer the equilibrium value of the following system

$$\begin{aligned} v(n) = 0.6 \cdot v(n-1) \\ p(n) = 0.13 \cdot v(n) + 0.87 \cdot p(n-1) + 25 \end{aligned}$$

with initial values $v(0) \approx 1441.67$, $p(0) = 3000$ and $p(1) = 3500$. However I've only seen equilibrium point calculations done on equations of the form

$$y(k+1)=ay(k)+b$$

which only has one function on the R.H.S. My system has two: $p(n-1)$ and $v(n-1)$.

I cannot find information regarding how to calculate the equilibrium of such multiple function system. How is this done?

Note:

My system can be formulated into the matrix form:

$$x(n) = A\cdot x(n-1) + f$$

Can this be used to find the equilibrium, eventhough $x(n)$ is a $2 \times 1$ vector?

If I calculate the equilibrium using this matrix form, then I get

$x(n) = A \cdot x(n-1) + f$
Let $x(n) = f(x(n-1))$
$f(x(n-1)) = A \cdot x(n-1) + f$
Let $f(x(n-1)) = x(n-1)^*$
$x(n-1)^* = A \cdot x(n-1)^* + f$
From which
$(A-I) \cdot x(n-1)^* + f = 0$
$x(n-1)^*=(A-I)^{-1} \cdot (-f)$

And Matlab gives:

>> I = eye(size(A))
>> inv(A-I)*(-f)

ans =

         0
  192.3077

which looks reasonable. And this also seems to follow equilibrium behaviour, since:

xstar = inv(A-I)*(-f)
for i = 1:99
    xstar = [xstar A*xstar(:,end)+f];
end

xstar =

  Columns 1 through 8

         0         0         0         0         0         0         0         0
  192.3077  192.3077  192.3077  192.3077  192.3077  192.3077  192.3077  192.3077

... and so on ...

---

Extra question:

What about inferring the stability of such matrix form?

Answer 1

You did a good job of translating the recurrence relation into matrix form. There is indeed a theory for stability of systems of linear recurrences, and it is critically tied to the eigenvalues of $A$. We have the system: $$x_n = Ax_{n-1} + f$$ with a fixed point that you've already calculated is $x^* = - (A-I)^{-1} f$.

One special case is when $1$ is an eigenvalue of $A$, in which case $A-I$ doesn't have an inverse and additionally the system has infinitely many fixed points (or possibly none). Let's ignore that case for the present, since it doesn't occur in your question.

To study stability, we just need to consider the difference $y_n = x_n - x^*$ between $x_n$ and the (unique) fixed point. A routine calculation shows that this results in the homogeneous recursion

$$y_n = Ay_{n-1} \implies y_n = A^n y_0.$$

So it boils down to understanding the behavior of powers of $A$, and naturally this is related to the eigenvalues of $A$. In the one-variable case (powers of a scalar $a$) this is very straightforward: if $|a|>1$ then the fixed point is a repeller, and if $|a|<1$ then it's an attractor. If $a=+1$ then we have many fixed points, and if $a=-1$ then it oscillates without getting closer to or farther from the fixed point.

In many way the analysis for the two-variable case is similar, but naturally a few more complications can arise. If $A$ is diagonalizable then by a change of basis we can take $A$ to be the diagonal matrix of eigenvalues $$A \sim \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix}.$$

The powers of this matrix are easily understood just as in the one-dimensional case: if $|\lambda_1|,|\lambda_2|<1$ then the fixed point is an attractor, and if $|\lambda_1|,|\lambda_2|>1$, then the point is a repeller. Of course it may also happen that $|\lambda_1|>1$ and $|\lambda_2|<1$: this case is called a saddle point, where most trajectories escape from the fixed point, but there is a basin of attraction (in this case a straight line) within which the point acts as an attractor.

Another interesting case that doesn't happen in one dimension is where $|\lambda_1| = |\lambda_2| = 1$ and $\lambda_1,\lambda_2$ are complex: this is similar to the $a=-1$ oscillation, except that here the orbit rotates around the fixed point in some sort of elliptical. Complex eigenvalues have this effect in general of creating rotation, though if $|\lambda_1| = |\lambda_2| < 1$ then this is more of a "spiralling inwards" (or outwards in the $>1$ case). Aside from this, the behaviour can be understood just as easily as the sequence $a^n$ for various choices of $a$.

Finally, it may happen that $A$ is not diagonalizable, in which case we have a repeated (real) eigenvalue $\lambda$ which by a change of basis corresponds to a Jordan block $$A \sim \begin{bmatrix}\lambda & 1 \\ 0 & \lambda\end{bmatrix}.$$

The behaviour of powers of $A$ in this form is in most cases the same as the one-dimensional case: if $|\lambda|>1$ then it's a repeller and if $|\lambda|<1$ it's an attractor. If $\lambda=1$ then the fixed point is not unique: there is an entire line of fixed points. For stability, if we start on this line then of course we stay on it, but if we start off of the line then it repels further away. Finally if $\lambda=-1$ then $$A^2 \sim \begin{bmatrix}1 & -2\\ 0 & 1\end{bmatrix} \sim \begin{bmatrix}1 & 1\\ 0 & 1\end{bmatrix}.$$ This case combines the oscillation of $a=1$ with the previous behavior: there is a line on which it oscillates, and off of that line it diverges away from that line, in an oscillating way like the sequence $(-1)^n n$.

Answer 2

This is not an answer but it is too long for a comment.

Making the problem more general $$v_n=a v_{n-1}\qquad p_n=b v_n+(1-b)p_{n-1}+c\qquad (0<a<1,0<b<1)$$ it is possible (even if tedious) to have closed form expressions for $v_n$ (this is the easy part) and for $p_n$ (this is the hard and tedious part).

Using limits, the funny result is, $$\lim_{n\to \infty } \, p_n=\frac c b$$