How to numerically solve difference equations in python

Question

I am trying to learn how to solve difference equations (also called recurrence relations) using python.

The problem in question is the equation

$x_{n+2} - 4x_{n+1} - x_{n} = 0$       where x_0 = 1 and x_1 = 1

Which outputs the sequence: n = 1, 1, 5, 21, 89, 377, ....

I have solved the problem by using math to find the general expression for this sequence, and then defining it as a function in python - and made it work (kind of).

However, my reason for posting this is that I think there are better ways to do this and that the solution I mentioned above is suboptimal.

After looking at a few similar examples, like how to numerically compute the Fibonacci sequence, i tried to generalize this method and modify it to the problem I'm working on hoping that it would work. And it kind of did, but not quite.

The solution I came up with is:

from numpy import zeros

N = 100
x = zeros(N+1, int)

x[0] = 1
x[1] = 1

for n in range(2, N+1):
    x[n] = x[n-2] + 4*x[n-1]
    print(f"x[{n}] = {x[n]}")

##produces wrong calculations after x[15]##

So my two questions are:

1. Is there a general and "solid" way of solving these kinds of problems? If so, does anyone have any examples they would like to share?

2. Why am I getting weird results after x[15]? And is there anyone who can help me fix it?

The out-print is supposed to be

Iteration    sequence output
===========================
x[0]        1
x[1]        1 
x[2]        5
x[3]        21
x[4]        89
x[5]        377
x[6]        1597
x[7]        6765 
…….
x[100]  137347080577163458295919672868222343249131054524487463986003968

And I’m getting:

x[0]        1
x[1]        1 
x[2]        5
x[3]        21
x[4]        89
x[5]        377
x[6]        1597
x[7]        6765
...
x[15] = 701408733
…….
x[16] = -1323752223
x[17] = -298632863
x[18] = 1776683621
x[19] = -1781832971
...
x[100] = -855830631

Answer 1

Numpy was developed for fast operations on arrays of numbers. This means that the type hint gets translated into a numpy-internal number type of fixed bit length. One can abuse numpy for operations on other object by defining dtype=object, but this will in general not be faster than list and iterator based solutions.

You define the type of the numpy zeros array as int, which apparently gets cast as 32 bit integers. Thus you get results in the limitations of the 32 bit integers with overflow folding into the negative range. Check with print(x.dtype), you should get int32 or, as with me in python3, int64.

Python's built-in integer objects that are used automatically in all integer operations have a multi-precision or bigint type, but that is not very compatible with numpy. Try

x=(N+1)*[0]

or list construction by appending as numpy-free alternative. Then with only this change your code produces

x[14] = 165580141
x[15] = 701408733
x[16] = 2971215073
x[17] = 12586269025
...
x[30] = 1779979416004714189
x[31] = 7540113804746346429
x[32] = 31940434634990099905
x[33] = 135301852344706746049
...
x[100] = 137347080577163115432025771710279131845700275212767467264610201