0

我正在尝试学习如何使用 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]##

所以我的两个问题是：

  1. 是否有解决此类问题的通用且“可靠”的方法？如果是这样，有没有人有任何他们想分享的例子？

  2. 为什么在 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
4

2 回答 2

2

差分方程只是递归关系。它们背后的数学可能非常棘手......找到伴随矩阵的基础......但你需要一个坚实的背景。对于那些东西，我建议你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

为了提高递归，您可以使用docs : memoizing callable 中的装饰lru_cache器，functools它最多可以保存最近调用的 maxsize。我感谢 Nachiket 提出这样的建议。

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

于 2021-09-16T08:29:16.327 回答
1

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

您将 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
于 2021-09-16T08:42:17.387 回答