Limit of a recursive relation from initial iterates.

Question

The query: How to find limits of a converging recursive sequence from some (you may take as many as needed) of the consecutive initial iterates?

Example: I have a recursive relation of the form $\vec{x}(n+3)=a_2 \vec{x}(n+2)+ a_1 \vec{x}(n+1)+ a_0 \vec{x}(n) +\vec{c}$, where $\vec{x}(k),\vec{c}\in \mathbb{R}^2$ and $a_2,a_1, a_0 \in \mathbb{R}$ are chosen such that the relation converges to a limit point, $\frac{\vec{c}}{1-(a_2+a_1+a_0)}$. If one knows the constants $a_2,a_1, a_0$ and $\vec{c}$, one could compute the limit point.

However, I am left with few (say 10) initial iterates of $\vec{x}(k)$. I need to find the limit point. Here is the approach I tried.

Step 1 : Estimate $a_2, a_1$ and $ a_0$ from the iterate values solving the following equation.

for k =1, 2 (or any two values)

$\vec{x}(k+5)-\vec{x}(k+4)$ = $a_2(\vec{x}(k+4)-\vec{x}(k+3))+a_1(\vec{x}(k+3)-\vec{x}(k+2))+a_0(\vec{x}(k+2)-\vec{x}(k+1))$

this is expected to give me four equations for three unknowns (I could generate as many equations as I need from the iterates varying $k$)

Here is the Problem: The system of equations is unsolvable (Matlab solver says the system is inconsistent.)

Step 2: Estimate $\vec{c}$ from the recursive relations.

Step 3: Estimate $x_{limit}$ using aforementioned problem.

Please find attached the sample Matlab code (if you wish to try out).

I observed that the condition number of the matrix formed for solving the linear equations mentioned above is large (which makes the system inconsistent, I guess). Why is it happening?

Any effort made in helping me solve this problem (or even identifying why this is happening) is highly appreciated.

Thanks in advance

%% Sample code


%% Calculating iterates using a sample recursive relation
a_2 = 0.1;
a_1 = 0.1;
a_0 = 0.3;

c=[0.2,1]';


x_limit = c/(1-(a_2+a_1+a_0));%-----(1)
% x_0 = 10*rand(size(c));
x_0 = zeros(size(c));
x_1 = a_0*x_0 +c;
x_2 = a_1*x_1+a_0*x_0 +c;

totla_iterates=200;
X_mat=zeros(length(c),totla_iterates);
X_mat(:,1) = x_0;
X_mat(:,2)=x_1;
X_mat(:,3)=x_2;
for i=4:totla_iterates
   X_mat(:,i)=a_2* X_mat(:,i-1)+a_1* X_mat(:,i-2)+a_0* X_mat(:,i-3)+c;%----(2)
end

%% Note that the recursive relation converges to the limit predicted by (1)

% Now use the first 10 iterates of (2) to predict the limit
%% Step 1: Estimate a_0,a_1
X_difmat=X_mat(:,2:end)-X_mat(:,1:end-1);

syms a2 a1 a0
i=1; % you may try varying the value of i or use more equations to find the soln.
% eqn1 =  X_difmat(1,i)*a2+X_difmat(1,i+1)*a1 +X_difmat(1,i+2)* a0 == X_difmat(1,i+3);
% eqn2 =  X_difmat(2,i)*a2+X_difmat(2,i+1)*a1 +X_difmat(2,i+2)* a0 == X_difmat(2,i+3);
% eqn3 =  X_difmat(1,i+1)*a2+X_difmat(1,i+2)*a1 +X_difmat(1,i+3)* a0 == X_difmat(1,i+4);
% [A,B] = equationsToMatrix([eqn1,eqn2,eqn3], [a2 a1 a0]);
eqn1 =  X_difmat(:,i)*a2+X_difmat(:,i+1)*a1 +X_difmat(:,i+2)* a0 == X_difmat(:,i+3);
eqn2 =  X_difmat(:,i+1)*a2+X_difmat(:,i+2)*a1 +X_difmat(:,i+3)* a0 == X_difmat(:,i+4);
[A,B] = equationsToMatrix([eqn1,eqn2], [a2 a1 a0]);

X=double(linsolve(A,B)); % note that I am unable to calculate a_1 and a_0 here
 disp(num2str(X)) % Ideally this should be X= a_2 a_1 and a_0 , which I am not getting.

```

Answer 1

This is an interesting problem which has a surprisingly simple solution. We start with a simple form of the problem. Suppose $\, x_{n+1} = a_0 x_n + c\,$ for all $\,n.\,$ We want to solve for $\,a_0,\,c\,$ and compute $\, L := c/(1-a_0).\,$ We solve the linear system $\, x_1 = a_0 x_0 + c,\, x_2 = a_0 x_1 + c\,$ and find that the expression for the limit $\,L_2\,$ (if it exists) is $$ L_2 = \frac{x_0 x_2 - x_1 x_1}{x_0 - 2x_1 + x_2}. $$ The numerator is the determinant of a $\,2 \times 2\,$ Hankel matrix formed using $\,(x_0, x_1, x_2).\,$ The denominator is the Total derivative of the numerator with all the partial derivatives replaced with $\,1.$ Notice that $\,L_2\,$ is exactly the result of Aitken's $\Delta^2$-quared process

This rational expression for $\,L_2\,$ naturally generalizes for linear recurrences with more terms. For example, suppose that $\, x_{n+2} = a_1 x_{n+1} + a_0 x_n + c\,$ for all $\,n\,$ and the limit $\, L := c/(1-a0-a1).\,$ Solving the linear system $\, x_2 = a_1 x_1 + a_0 x_0 + c,\, x_3 = a_1 x_2 + a_0 x_1 + c,\, x_4 = a_1 x_3 + a_0 x_2 + c\,$ gives the expression for the limit $\,L_3\,$ (if it exists) as $$ L_3 = \frac{ x_0 x_2 x_4 + 2 x_1 x_2 x_3 - x_2^3 - x_0 x_3^2 - x_1^2 x_4 } { (x_0 - 2 x_1 + x_2) (x_2 - 2 x_3 + x_4) - (x_1 - 2 x_2 + x_3)^2}. $$

The reason for this general result is that the numerator of the limit $\,L\,$ is $\,c\,$ and $\,c=0\,$ is equivalent to the Hankel determinant of the homogeneous linear system being zero. If the denominator $\,(1 - a_0 - ... - a_k) = 0,\,$ then the 2nd difference of the $\,\{x\}\,$ sequence satisfies a linear homogeneous system of equations and is equivalent to a Hankel determinant being zero.

Of course, this assumes that the limit exists and exact calculations are used. One issue is that the denominator and/or the numerator could be zero. For example, if $\, x_{n+1} = x_n\,$ then $\, L_2 = 0/0 \,$ which gives no information about the limit of the constant sequence. If $\,\{x\}\,$ is a geometric sequence given by $\, x_{n+1} = a_0 x_n\,$ then $\, L_2 = 0/(x_0(1-a_0)^2)\,$ which implies $\, L_2 = 0\,$ if $\,x_0 \ne 0\,$ and $\, a_0 \ne 1,\,$ but the limit is zero only if $\, |a_0|<1.\,$ The other issue is loss of significance in doing inexact arithmetic.

For testing purpose, I wrote the following PARI/GP code:

{doit(n, m=1) = for(k=0, n+m-1, x[n+k+1] = c + sum(i=1,n, x[k+i]*a[i]))};
{L3(k=1) = my(x0 = x[k+0], x1 = x[k+1], x2 = x[k+2], x3 = x[k+3], x4 = x[k+4],
   y0 = x0 - 2*x1 + x2, y1 = x1 - 2*x2 + x3, y2 = x2 - 2*x3 + x4);
   (x0*x2*x4 + 2*x1*x2*x3 - x2^3 - x0*x3^2 - x1^2*x4)/(y0*y2 - y1*y1)};
default(realprecision, 9);
ni = 20; x = vector(ni+4); a = [0.2,0.1]; c = 0.02; x[1] = 0.5; x[2] = 0.3;
   doit(2, ni); print(x); print(vector(ni, k, L3(k)));

The resulting output is:

[0.500000000, 0.300000000, 0.150000000, 0.0950000000, 0.0595000000, 0.0449500000, 0.0363950000, 0.0326295000, 0.0305419500, 0.0295800950, 0.0290663995, 0.0288226589, 0.0286955458, 0.0286340864, 0.0286025178, 0.0285870690, 0.0285792105, 0.0285753349, 0.0285733756, 0.0285724045, 0.0285719156, 0.0285716725, 0.0285715504, 0.0285714895]
[0.0285714280, 0.0285714283, 0.0285714283, 0.0285714288, 0.0285714246, 0.0285714317, 0.0285714052, 0.0285716963, 0.0285716927, 0.0285761767, 0.0285658400, 0.0286175962, 0.0286606325, 0.0258589033, 0.0269360081, 0.E-1, 0.E0, -1.34551706, -10.0552618, 33.4275327]

You can see the original sequence converging, but the approximations to the limit appear to diverge. This is the result of using only $9$ digits of precision. If the number of digits of precision is increased to $19$, the problem goes away.

Answer 2

If your iterates are indeed coming from the given recurrence, the system must be compatible. In fact, there are four unknown constants and you should work with four equations. Then the coefficients so obtained should let you retrieve the next iterates exactly.

If your problem is numerically unstable, you might be facing significant numerical errors, and it can make sense to work with an overdetermined system and solve it in the least-squares sense (or similar).

If the recurrence is only hypothesized or accepted as an approximation or if the data is noisy, then perforce the systems will be incompatible and you need to resort to a best-fit.

---

By the way, there is no reason to treat the independent term $\vec c$ differently than the other unknown parameters.