如何数值求解 python 中的差分方程

Question

我正在尝试学习如何使用 python 求解差分方程（也称为递归关系）。

问题是方程

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

输出序列：n = 1, 1, 5, 21, 89, 377, ....

我已经解决了这个问题，方法是使用数学找到这个序列的一般表达式，然后将它定义为 python 中的一个函数 - 并使其工作（有点）。

但是，我发布这个的原因是我认为有更好的方法可以做到这一点，而且我上面提到的解决方案不是最优的。

看了几个类似的例子，比如如何用数字计算斐波那契数列，我试图推广这种方法并将其修改为我正在处理的问题，希望它能奏效。它有点像，但不完全是。

我想到的解决办法是：

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]##

所以我的两个问题是：

是否有解决此类问题的通用且“可靠”的方法？如果是这样，有没有人想分享任何例子？
为什么我在 x[15] 之后得到奇怪的结果？有没有人可以帮我解决这个问题？

外刊应该是

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

我得到：

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

差分方程只是递归关系。它们背后的数学可能非常棘手……找到伴随矩阵的基础，……但你需要扎实的背景知识。对于那些东西，我建议你sympy这是一个用于符号操作的数学包。

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))

输出

1
5
21
89
377
1597
6765
28657
121393
514229
2178309
9227465
39088169
165580141
701408733
2971215073
12586269025
53316291173
225851433717

要提高递归，您可以使用 functools、docs 中的装饰器 lru_cache：memoizing 可调用函数最多可保存最近调用的最大大小。感谢 Nachiket 的建议。

关于你的第二个问题：你忘记添加错误日志：RuntimeWarning: overflow encountered in long_scalars。如果您使用 np.int64，您可能达到的最大整数是 9223372036854775807。

Answer 2

Numpy 是为对数字数组进行快速运算而开发的。这意味着类型提示被转换为固定位长度的 numpy 内部数字类型。可以通过定义 dtype=object 滥用 numpy 对其他对象进行操作，但这通常不会比基于列表和迭代器的解决方案更快。

您将 numpy zeros 数组的类型定义为 int，它显然被转换为 32 位整数。因此，您会得到 32 位整数限制的结果，溢出折叠到负数范围内。检查 print(x.dtype)，你应该得到 int32 或者像我一样 python3，int64.

Python在所有整数运算中自动使用的内置整数对象具有多精度或bigint类型，但与numpy不太兼容。尝试

x=(N+1)*[0]

或通过附加作为无 numpy 的替代列表构造。然后仅通过此更改，您的代码就会生成

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

如何数值求解 python 中的差分方程

How to numerically solve difference equations in python

python

math

numpy

numerical-methods

difference-equations