Is the iterative process convergent? My task is to investigate the convergence of the following iterative process:
$$\frac{\mathrm{x}^{n+1}-\mathrm{x}^{n}}{\tau} + \frac{1}{4}A(3\mathrm{x}^{n+1}+2\mathrm{x}^{n}-\mathrm{x}^{n-1}) = \mathrm{b}$$
I get an expression for $\mathrm{x}^{n+1}$:
$$ \mathrm{x}^{n+1}= \big(\mathrm{b} - \big(\frac{1}{2}A - \frac{1}{\tau}\big)\mathrm{x}^{n} + \big(\frac{1}{4}A \big)\mathrm{x}^{n-1}\big) \big(\frac{1}{\tau}E + \frac{3}{4}A\big)^{-1} $$
What to do next is not clear to me, in particular, how to look for a transition matrix
$B$? How to get rid of $\mathrm{x}^{n-1}$?
I think we can make a tricky replacement here, but I can't think of one.
Сan anyone help with the solution?
 A: This is an implicit multi-step method for the differential equation $\dot x+Ax=b$. If $A$ is positive (real part of all eigenvalues positive) and the step size inside the stability region (scaled by $-A$), the iteration should converge to $x=A^{−1}b$.
To explore that stability region, use the exponential trial solution $x(t)=ce^{-λt}$ for the homogeneous equation, so that $x^n=ce^{-λnτ}=cq^n$. For the exact solution one knows that $λ=A$ or an eigenvalue of $A$.
As the homogeneous equation is homogeneous, set $c=1$ and replace $x^n=q^n$, left index, right power, divide by the lowest power of $q$ and find the two solutions of the resulting quadratic equation in $q$. These solutions will depend on $τ$. For convergence you need both roots inside the unit circle.
import sympy
q, A, tau = sympy.symbols("q A \\tau")
eq = sympy.Eq(0, (q**2-q)/tau+A/4*(3*q**2+2*q-1))
rts = sympy.solve(eq,q)
print(sympy.latex(rts))

$$
\left [ \frac{- A \tau - 2 \sqrt{A^{2} \tau^{2} + 1} + 2}{3 A \tau + 4}, \quad \frac{- A \tau + 2 \sqrt{A^{2} \tau^{2} + 1} + 2}{3 A \tau + 4}\right ]
$$
This does not appear to be immediately helpful.
Another possibility is to trace the unit circle with $q$ and plot the contours that the value $z=\tau A$ (if scalar, else $z=λτ$) can take
q = np.exp(1j*np.pi*np.linspace(-0.75,0.75,300))
R = lambda q: -4*q*(q-1) / (3*q-1)/(q+1)

z = R(q)
plt.plot(z.real, z.imag, 'b')
for k in range(1,4):
    z = R(0.995**k*q)
    plt.plot(z.real, z.imag, '.g', ms=2)
plt.grid(); plt.show()

The blue line is the boundary, the green dots indicate the interior of the stability region.

This at least gives the idea that stability is to be had independent of $\tau>0$ for all $A>0$, and if $A$ is a matrix, if all eigenvalues have a positive real part.
A: The most straightforward approach for this type of convergence analysis should be as follows

Let $(x^n,n \in \mathbb{N})$  be a sequence of vector in $\mathbb{R}^{d \times 1}$, $B$ and $C$ be two matrices in $\mathbb{R}^{d \times d}$, $b$ be any vector in $\mathbb{R}^{d \times 1}$.
Suppose that  the following recursive relation holds for all $n$
$$x^{n+1}=Bx^n+Cx^{n-1}+b$$
We define $$v^n=
  \begin{bmatrix}
    x^n \\
     x^{n-1}
  \end{bmatrix} \in \mathbb{R}^{ 2d \times 1}$$ and 
$$T=  \begin{bmatrix}
    B & C\\
    \text{Id}_n & 0
  \end{bmatrix} \in \mathbb{R}^{2d \times 2d} \quad \tilde{b}= 
  \begin{bmatrix}
    b \\
     0
  \end{bmatrix} \in \mathbb{R}^{ 2d \times 1}$$
then
$$v^{n+1}=Tv^n+\tilde{b}$$
This type of recursive sequence converges for all possible initial values $(v^1,\tilde{b})$ iff the absolute of every eigenvalue of $T$ is smaller than $1$

Hence if you can show the later in the previous statement, you can show your desired convergence.
