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:
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?
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
CodePudding user response：
Difference equations are just recursive relationships. The mathematics behinds them can be quite tricky... finding the basis of the companion matrix, ... but you need a solid background. For that stuffs I suggest you sympywhich is a mathematics package for symbolic manipulation.
import functools
@functools.lru_cache
def diff_eq(n):
if n == 0 or n == 1: return 1
return 4*diff_eq(n-1) diff_eq(n-2)
for i in range(20):
print(diff_eq(i))
Output
1
5
21
89
377
1597
6765
28657
121393
514229
2178309
9227465
39088169
165580141
701408733
2971215073
12586269025
53316291173
225851433717
To boost the recursion you can use the decorator lru_cache from the functools, docs: memoizing callable that saves up to the maxsize most recent calls. I am thankful to Nachiket for such suggestion.
For your second question: you forgot to add the error log: RuntimeWarning: overflow encountered in long_scalars. If you use np.int64 the maximum integer that you may achieve is 9223372036854775807.
CodePudding user response：
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