如何在Python中数值求解差分方程（递推关系）及异常结果修复

问题2：代码异常的原因与修复方案 #

你遇到的计算错误是numpy默认整型溢出导致的：
numpy的int类型默认是平台相关的32位有符号整型，可存储的最大值为2^31 -1 = 2147483647，计算x[16]时结果为4*701408733 + 165580141 = 2971215073，已经超过32位有符号整型的上限，触发整数溢出后高位符号位被改写，就会出现负数结果，后续计算也会全部出错。

修复方法非常简单，直接用Python原生列表存储序列即可，Python的内置int支持任意精度，不会出现溢出问题：

N = 100
x = [0]*(N+1)
x[0] = 1
x[1] = 1

for n in range(2, N+1):
    x[n] = 4 * x[n-1] + x[n-2]

# 打印结果
print("迭代次数    序列输出值")
print("="*27)
for n in range(N+1):
    print(f"x[{n:<6d} {x[n]}")

如果一定要用numpy，可将数组的dtype设为object，用Python原生整型存储：

from numpy import zeros
N = 100
x = zeros(N+1, dtype=object)
x[0] = 1
x[1] = 1
for n in range(2, N+1):
    x[n] = 4 * x[n-1] + x[n-2]

问题1：线性递推关系的通用求解方案 #

针对这类常系数线性齐次递推关系，有三类常用的可靠实现方案，适用不同场景：

1. 逐项迭代法 #

就是你当前使用的思路，实现最简单，适合需要输出前N项所有值的场景，时间复杂度O(N)，空间复杂度可以优化到O(1)（只存前两项即可，不需要存整个序列）：

# 空间优化版迭代，仅需存储前两项
def get_seq(N):
    if N == 0:
        return 1
    a, b = 1, 1
    print("x[0] = 1")
    print("x[1] = 1")
    for n in range(2, N+1):
        c = 4*b + a
        print(f"x[{n}] = {c}")
        a, b = b, c
    return b

2. 矩阵快速幂法 #

如果需要计算单个非常大的n对应的项（比如n=1e18），迭代法效率太低，可以用矩阵快速幂，时间复杂度仅为O(logN)。二阶递推x_{n+2} = p*x_{n+1} + q*x_n的通用矩阵形式为：
$$
\begin{bmatrix} x_{n+1} \ x_n \end{bmatrix} = \begin{bmatrix} p & q \ 1 & 0 \end{bmatrix}^n \cdot \begin{bmatrix} x_1 \ x_0 \end{bmatrix}
$$
对应你的问题p=4，q=1，实现示例如下：

def matrix_mult(a, b):
    # 2*2矩阵乘法
    return [
        [a[0][0]*b[0][0] + a[0][1]*b[1][0], a[0][0]*b[0][1] + a[0][1]*b[1][1]],
        [a[1][0]*b[0][0] + a[1][1]*b[1][0], a[1][0]*b[0][1] + a[1][1]*b[1][1]]
    ]

def matrix_pow(mat, power):
    # 矩阵快速幂
    result = [[1,0], [0,1]] # 单位矩阵
    while power > 0:
        if power %2 ==1:
            result = matrix_mult(result, mat)
        mat = matrix_mult(mat, mat)
        power = power //2
    return result

def get_seq_item(n):
    if n ==0:
        return 1
    if n ==1:
        return 1
    mat = [[4,1], [1,0]]
    mat_p = matrix_pow(mat, n-1)
    return mat_p[0][0] *1 + mat_p[0][1] *1

3. 通项公式法 #

你已经推导过闭式解，直接代入计算即可，仅适合n较小的场景，n过大会出现浮点数精度丢失问题。

内容的提问来源于stack exchange，提问作者fredericoamigo