Recurrence relation is not numerically stable.

Question

I have to solve the following Recurrence relation both analytically and numerically. Solving it analytically gives the answer $2a(\frac{1}{2})^n$ which is stable for every value of a, but when I try to evaluate the relation in numerically in python, it gives very large numbers for big $n$ I know it has something to do with the value of a. When I try to calculate the relation with a = 1 it is stable, but calculating it with $a = \frac{1}{3}$ gives very strange results. $$X_n = \frac{5}{2}X_{n-1} - X_{n-2}$$ $$X_1 = a, \quad X_2 = \frac{a}{2}$$I don't understand what causes the error propagation.

My code:

n = []
n.append(1/3)
n.append(1/6)
answer = 0
for i in range(2,1000,1):
    n.append(5/2 *n[i - 1] - n[i - 2])
print(n[999])

Answer 1

Starting from @Ian' comment, consider $$X_n = \frac{5}{2}X_{n-1} - X_{n-2}\qquad \text{with}\qquad X_1 = a\quad \text{and}\quad X_2 = \frac{a}{2}+\color{red}{\epsilon}$$ The solution is $$X_n=\frac a {2^{n-1}}+\color{red}{\frac{2^n }{3}\left(1-\frac 1{2^{2(n-1)}} \right)\epsilon }$$

If you make $n=1000$ and $\epsilon=10^{-100}$, the red term is $3.57\times 10^{200}$

Answer 2

The general solution is of the form $X_n=c_1 2^n + c_2 2^{-n}$. The analytical solution has $c_1=0$ no matter what $a$ is. But after one step of the numerical recurrence, the result, call it $\tilde{X}_3$, may not be exactly $a/4$, because of roundoff. If this happened and then you reinitialized the analytical recurrence with $\tilde{X}_3$ and $X_2$ as the initial data then there would be a small $c_1$, whose contribution would grow exponentially. It should therefore not be a huge surprise that the actual numerical trajectory does something similar.