# Bound on matrix product $\begin{bmatrix} 1+\frac{1}{n} & -1 \\ 1 & 0 \end{bmatrix}\cdots\begin{bmatrix} 1+\frac{1}{2} & -1 \\ 1 & 0 \end{bmatrix}$

Let matrices $$A_n = \begin{bmatrix} 1+\frac{1}{n} & -1 \\ 1 & 0 \end{bmatrix}$$ and the matrix product $$X_n = A_n A_{n-1} \cdots A_3 A_2 \quad .$$ Determine a bound on $$|(X_n)_{1,1}|$$, independent of $$n$$.

The actual form of the product gets pretty involved. It appears from numerical simulations that $$|(X_n)_{1,1}| \le |(X_2)_{1,1}| = 1.5$$ but I wouldn't know how to prove it. A less tight bound, independent of $$n$$, would be fine too.

• have you tried playing around with matrix decomposition? maybe there's some easy way to bound $S_{n-1} \cdot S_{n}^{-1}$, where $A_n = S_n^{-1} \Lambda_n S_n$. If so, possibly we could just take this bound (hopefully uniform) and then use the product of first eigenvalues. Commented Apr 5, 2020 at 16:57
• The continuous analogue for the evolution of the components of $X$ (the recurrence relation) would be that of a harmonic oscillator $y'' + w^2y = 0$ with a spring constant that changes with time (it gets stiffer) giving rise to a vibrartion frequency $w^2 = (1- 1/(t+1))$. For the continuous problem one can use energy methods to show that the energy of this system is bounded which implies a bound on $X$. The initial velocity is $0$ and the initial excursion is $y = 1.5$ so this is the bound I would expect to hold. Commented Apr 5, 2020 at 18:30
• @Winther Do you know of a way to also solve this (i.e. find a bound) by perturbation of the solution with fixed $\omega$? Commented Apr 5, 2020 at 18:38
• @Winther Thanks. I searched around for energy method and found a treatment at math.stackexchange.com/questions/3604480 . Commented Apr 5, 2020 at 18:38
• Yes the same method can be applied to the discrete case, its just a bit more messy as you can see in that linked question (which luckily is exactly the same problem as this one reduces to - even the initial conditions match). Commented Apr 5, 2020 at 18:45

TLDR; We can get $$x_n^2 \leq \frac{8}{3}$$ directly and with some computation we can establish the best bound $$x_n^2 \leq \frac{9}{4}$$ for all $$n$$. The argument below can be made much shorter, I just presented it as it is as to show the thought process behind it.

With $$X_n = \pmatrix{a_n & b_n \\ c_n & d_n}$$ then both $$a_n$$ and $$b_n$$ satisfy the recurrence $$x_{n+1} = \left(1 + \frac{1}{n}\right)x_n - x_{n-1}$$ with $$c_n = a_{n-1}$$ and $$d_n = b_{n-1}$$ (and are linked via the determinant constraint $$a_nb_{n-1} - b_na_{n-1} = 1$$).

We can write this recurrence on the more suggestive form $$\frac{x_{n+1}+x_{n-1}-2x_n}{h^2} + \left(1 - \frac{1}{t_n}\right)x_n = 0$$ where the step-size $$h=1$$ and $$t_n = hn$$. This can also be written in terms of the difference operator $$\Delta x_n \equiv \frac{x_{n+1}-x_n}{h}$$ as $$\Delta^2 x_{n-1} + w_n^2x_n = 0$$ where $$w_n^2 = 1 - \frac{1}{t_n}$$. From this form its easily seen to be a discretization of the harmonic oscillator ODE $$\frac{d^2x(t)}{dt^2} + w^2(t) x(t) = 0$$ with a time-dependent frequency $$w^2(t) = 1-\frac{1}{t}$$. The other way we can see the relation with an ODE is to notice that the product of the matrices you have is what you would naturally get when you try to evolve a system of two coupled ODEs with an explicit method.

The asymptotic solution is just $$x(t) = A \cos(wt + \phi)$$ just as the asymptotic solution to the discrete equation (ignoring the $$1/n$$ term) is $$\cos(2\pi n / 6 + \phi)$$.

The analogous problem of bounding the components of $$X_n$$ is in the continuous setting the same as bounding $$x(t)$$. To do so the standard method is to consider the evolution of the energy of the system. When $$w(t)$$ is a constant then we have that the energy $$\frac{1}{2}\dot{x}^2 + \frac{1}{2}w^2x^2$$ is constant. Defining the energy to be $$E(t) = \frac{1}{2}\dot{x}^2 + \frac{1}{2}w^2(t)x^2$$ then we find that $$\dot{E} = w\dot{w}x^2$$ and as long as $$\dot{w} \geq 0$$ and $$w$$ is bounded then $$E$$ is increasing and $$\dot{E} \leq \frac{2\dot{w}(t)}{w}E \implies 1 \leq \frac{E(t)}{E(t_0)} \leq \frac{w^2(t)}{w^2(t_0)}$$ so $$E$$ is bounded and therefore $$x$$ is bounded by $$\frac{\sqrt{2E_{\rm max}}}{w}$$.

The discrete problem is meant to approximate the continuous one so we would expect the same type of method to work (but its discrete so it will not be as clean as above). We will try to define an energy similar to the one above. The natural choice for the kinetic energy is

$$E^K_n = \frac{1}{2}(\Delta x_{n-1})^2$$

which satisfy

$$\Delta E^K_n = \frac{x_{n+1} - x_{n-1}}{2h}\cdot \frac{x_{n+1} + x_{n-1} - 2x_n}{h^2}\\ = (\Delta^{\rm symm} x_n) \cdot (\Delta^2 x_{n-1}) = -(\Delta^{\rm symm} x_n)w_n^2x_n$$ where $$\Delta^{\rm symm} x_n$$ is the symmetric difference formula.

For the potential energy we must do a slight twist and define it as $$E^P_n = \frac{1}{2}F_n\left(\frac{x_n+x_{n-1}}{2}\right)^2$$ where $$F_n = \frac{w_n^2}{1-w_n^2h^2/4}$$. This leads to $$\Delta E^P_n = \frac{1}{2}\frac{F_{n+1}-F_n}{h}\left(\frac{x_n+x_{n+1}}{2}\right)^2 + (\Delta^{\rm symm} x_n)w_n^2 x_n$$ The particular form for $$F_n$$ we choose (instead of the more natural $$F_n = w_n^2$$) is needed to get the last term on the desired form as to cancel with the term from the kinetic energy (but $$F_n \to w_n^2$$ as $$h\to 0$$ so its indeed consistent with the continuous definition). The total energy satisfies $$\Delta E_n = \frac{1}{2}(\Delta F_n)\left(\frac{x_n+x_{n-1}}{2}\right)^2$$ and if $$F_n$$ is an increasing function (which it is here) then $$\Delta E_n \geq 0$$ and $$\Delta E_n \leq \frac{\Delta F_n}{w_n^2}E_n$$ We can now do the analogous thing as in the continuous case (i.e. apply Gronwalls lemma). This gives us $$E_M \leq E_N \leq E_M\exp\left[\sum_{n=M}^{N-1}\frac{\Delta F_n}{w_n^2}h\right]$$ Asymptotically we have $$E_M \leq E_N \lesssim E_M\exp\left[\sum_{n=M}^{N-1}\frac{16}{9n^2}\right] \sim E_M\exp\left[\frac{16}{9M}\right]$$ so $$E_n$$ converges to $$E_*$$ and we can bound it accurately by summing to high enough $$M$$ (and computing $$E_M$$ numerically).

However, we can find a better bound by noting that in the continuous case the final inequality says that $$E/w^2$$ is decreasing. This motivates considering the quantity $$Z_n = \frac{E_n}{F_n}$$ which we find satisfy

$$\Delta Z_n = -\frac{\Delta F_n E_n^K}{F_nF_{n+1}} \leq 0$$ which implies $$E_n \leq F_n\frac{E_M}{F_M}\leq F_\infty\frac{E_M}{F_M}$$ where $$F_\infty = \lim_{n\to\infty}F_n = \frac{4}{3}$$.

To complete this we need to bound $$x_n$$ in terms of $$E_n$$. The total energy $$2E_n = (x_n - x_{n-1})^2 + \frac{F_n}{4}(x_n + x_{n-1})^2$$ is just an ellipse in phase space so if we take $$x_n - x_{n-1} = \sqrt{2E_n}\cos(\phi_n)$$ and $$x_n + x_{n-1} = \sqrt{2E_n}\frac{2}{\sqrt{F_n}}\sin(\phi_n)$$ for some angle $$\phi_n$$ it follows that $$2x_n = \sqrt{2E_n}(\cos(\phi_n) + \frac{2}{\sqrt{F_n}}\sin(\phi_n))$$ and computing the maximum of the right hand side we get the same type of bound as in the continuous case

$$|x_n| \leq |x_n|^{\rm Max} = \frac{\sqrt{2E_n}}{w_n}$$

and by using the bound above

$$x_n^2 \leq \frac{2E_M}{F_M(1-w_\infty^2/4)} = \frac{8}{3}\frac{E_M}{F_M}$$

where $$M$$ is any integer. Since $$E_1/F_1 = 1$$ we obtain $$x_n^2 \leq \frac{8}{3}$$ directly. Computing the first few values of $$x_n$$ numerically we can establish $$|x_n| \leq \frac{3}{2}$$ and since $$x_2 = \frac{3}{2}$$ we get $$|x_n| \leq \frac{3}{2}$$ for all $$n$$. It might be possible to get a slightly better asymptotical bound is needed.

The bounds above (apart from the numbers quoted) is general and applies to similar recurrences $$x_{n+1} = (2-w_n^2h^2)x_n - x_{n-1}$$ corresponding to products of matrices $$A_n = \pmatrix{2-w_n^2h^2 & -1 \\ 1 & 0}$$ as long as $$w_n$$ is increasing and converges as $$n\to\infty$$.

• The answer is so nice. Commented Apr 9, 2020 at 6:59
• amazing. It even kind of explains where the special number $6$ was coming from in the simulation by @rogerl Commented Apr 9, 2020 at 19:44

I hope I'm not wrong, but since no one has said this yet, it seems too good to be true.

What I did was just write a recurrence relation between the coefficients: $$\begin{cases} (X_{n+1})_{1,1}=(1+\frac{1}{n+1})(X_{n})_{1,1}-(X_{n})_{2,1} \\ (X_{n+1})_{2,1}=(X_{n})_{1,1} \end{cases}$$

You then get a $$2$$nd order recurrence relation with non-constant coefficients, where you can factor the operator and solve explicitly. Am I looking at this wrong? I'll delete my answer if so.

• Beat me to it (+1). I think that is correct. Unfortunately the resulting system has a pretty complicated solution, at least according to Mathematica. Commented Apr 5, 2020 at 17:16
• Good point: you can write $(X_{n+1})_{1,1}=(1+\frac{1}{n+1})(X_{n})_{1,1}-(X_{n-1})_{1,1}$ or simpler write $c_{n+1}=(1+\frac{1}{n+1})c_{n}- c_{n-1}$, yet I don't see how you would factor it? Commented Apr 5, 2020 at 17:28

This is not an answer, but here is a graph of $$(X_n)_{1,1}$$ for $$n$$ up to 30000. This is a plot of the points; they only appear to be lines. Each apparent line represents a plot of $$(X_n)_{1,1}$$ for a particular congruence class of $$n$$ modulo $$6$$. Top to bottom at $$30000$$, the congruence classes are $$0$$, $$1$$, $$5$$, $$2$$, $$4$$, $$3$$.

Your conjecture seems plausible; the behavior of the sequence is pretty interesting.

• There are a number of graphs - you said you wanted to display $(X_n)_{1,1}$ ? Commented Apr 5, 2020 at 16:53
• There only appear to be six branches in the sequence. Each of these apparent lines corresponds to a sequence $(X_n)_{1,1}$ for a particular value of $n$ modulo 6. Commented Apr 5, 2020 at 16:55
• Intersting observation. Without the term $1/n$, you would also get 6 branches which would then be constants. Maybe it's possible to compute the analytic behavior of the branches separately. Commented Apr 5, 2020 at 17:00

Computing the eigenvalues from the first 4000 products we have the graphics

In Red the eigenvalues location into the complex plane

Follows the script to obtain it

Pc[list_] := Module[{},
x1 = Re[list[[1]]];
y1 = Im[list[[1]]];
x2 = Re[list[[2]]];
y2 = Im[list[[2]]];
Return[Show[Graphics[{Red, Point[{x1, y1}]}], Graphics[{Red, Point[{x1, y1}]}]]]
]

A[n_] := N[{{1 + 1/n, -1}, {1, 0}}]
A0 = A[2];
n = 4000;
list = {Pc[Eigenvectors[A0]]};
For[i = 3, i <= n, i++,
A0 = A0.A[i];
AppendTo[list, Pc[Eigenvectors[A0]]]
]

Show[Join[list, {Graphics[Circle[]]}], Axes -> True, PlotRange -> All]