如何在Python中解决递归关系

Question

我正在尝试写代码来给一个递归关系算出数字答案。这个关系其实很简单，定义如下。变量 x 是一个整数。

• 当 i 大于 0 且小于 x 时，p(i) = p(i+2)/2 + p(i-1)/2
• p(0) = p(2)/2
• 当 i 大于等于 x 时，p(i) = 1

这段代码里也有这个内容。

from __future__ import division
def p(i):
    if (i == 0):
        return p(2)/2
    if (i >= x):
        return 1
    return p(i-1)/2+p(i+2)/2


x = 4
#We would like to print p(0)  for example.

不过，这样写并不能真正计算出 p(0)。那在 Python 中该怎么做呢？

---

有没有可能建立一个同时方程的系统，然后用 numpy.linalg.solve 来解决它呢？

Answer 1

我有点困惑，因为你的代码看起来应该正是这样做的。

def p(i):
    x = 4 # your constant should be defined in-function

    if (i == 0):
        return p(2)/2
    elif (i >= x):
        return 1
    return p(i-1)/2+p(i+2)/2

这里最大的问题在于你的递归。对于 p(1) 来说，它会这样做：

p(0)/2                    +                    p(3)/2
p(2)/2                    +           p(2)/2 + p(4)/2
p(1)/2                    +              p(1)/2 + 1/2
# each side of the operator is now the same as the original problem!
# that's a sure sign of infinite recursion.

你期望输出是什么呢？

Answer 2

这里的问题是，无论你从哪里开始，都会陷入无限递归，因为递归并不是很明显，而是变成了需要解决的线性方程组。如果这是一个需要用Python来解决的问题，我会用Python来计算这个方程组的系数，然后用克拉默法则来解它。

补充一下：你的未知数是 p(0), ..., p(x-1)。一开始就可以得到一个系数行向量是 (1, 0, -1/2, 0, ..., 0)（因为 p(0)-p(2)/2=0），其他的都是类似 (..., -1/2, 1, 0, -1/2, ... ) 的形式。总共有 x-1 个这样的方程（对应 p(1), ..., p(x-1)），所以这个方程组要么有唯一解，要么没有解。直观上看，似乎总是应该有一个唯一解。

最后两个方程会是唯一的，因为它们会涉及到 p(x) 和 p(x+1)，所以这些项会被省略；然后克拉默法则右侧的列向量应该是 (0, 0, ..., 0, 1/2, 1/2)，我认为是这样的。

Numpy 支持矩阵运算。

Answer 3

这不是对提问的直接回答，但这个页面在谷歌搜索“用Python解决递推关系”时排名第一，所以我决定写个回答。

如果你有一个线性递推关系，并且想找到递推公式，可以使用Sympy库里的 find_linear_recurrence 函数。举个例子，假设你有以下这个数列：0, 1, 3, 10, 33, 109, 360, 1189, 3927, 12970。然后下面的代码就可以生成这个递推关系：

import sympy
from sympy.abc import n
L = [0, 1, 3, 10, 33, 109, 360, 1189, 3927, 12970]
print(sympy.sequence(L, (n, 1, len(L))).find_linear_recurrence(len(L)))

输出结果是：

[3, 1]

所以你可以知道 A(n) = 3*A(n-1) + A(n-2)。

Answer 4

你说得对，这个问题可以用线性代数来解决。下面我做了一个简单的硬编码翻译。你的方程从 p(0) 到 p(3) 被重新排列，使得右边等于 0。而 p(4) 和 p(5) 作为递推关系中的基础情况，右边则是 =1。

• -p(0) + p(2)/2 = 0

• p(i-1)/2 - p(i) + p(i+2)/2 = 0 适用于 i 大于 0 且小于 x

• p(i) = 1 如果 i 大于等于 x

这里是为 n=4 硬编码的程序

import numpy
a=numpy.array([[-1,   0, 0.5,  0,   0,   0], # 0
               [0.5, -1,   0,0.5,   0,   0], # 1
               [0,  0.5,  -1,  0, 0.5,   0], # 2
               [0,    0, 0.5, -1,   0, 0.5], # 3
               [0,    0,   0,  0,   1,   0], # 4
               [0,    0,   0,  0,   0,   1], # 5
              ])
b=numpy.array([0,0,0,0,1,1])
# solve ax=b
x = numpy.linalg.solve(a, b)
print x

编辑，这里是程序化构建矩阵的代码，只测试了 n=4！

n = 4

# construct a
diag = [-1]*n + [1]*2
lowdiag = [0.5]*(n-1) + [0]*2
updiag = [0.5]*n
a=numpy.diag(diag) + numpy.diag(lowdiag, -1) + numpy.diag(updiag, 2)

# solve ax=b
b=numpy.array([0]*n + [1]*2)
x = numpy.linalg.solve(a, b)

print a
print x[:n]

这个输出结果是

[[-1.   0.   0.5  0.   0.   0. ]
 [ 0.5 -1.   0.   0.5  0.   0. ]
 [ 0.   0.5 -1.   0.   0.5  0. ]
 [ 0.   0.   0.5 -1.   0.   0.5]
 [ 0.   0.   0.   0.   1.   0. ]
 [ 0.   0.   0.   0.   0.   1. ]]
[ 0.41666667  0.66666667  0.83333333  0.91666667]

与您在问题下评论中的解决方案相匹配。