Sympy rsolve函数返回错误结果：求解递推关系仅得到常数而非含n的表达式

问题背景 #

我最近在求解一个递推关系时遇到了麻烦，递推式是这样的：
$$a_n = 0.996 \cdot a_{n–1} + 0.004 \cdot a_{n–10} + 0.04$$
初始条件很简单：所有 $n ≤ 9$ 时，$a_n = 0$。

我用Sympy的rsolve函数尝试符号求解，写了这段代码：

import sympy as sp

n = sp.symbols('n')
y = sp.Function('y')
recurrence = sp.rsolve(
    y(n) - sp.Rational(0.996) * y(n-1) - sp.Rational(0.004) * y(n-10) - sp.Rational(0.04),
    y(n), 
    {
       y(0): 0,
       y(1): 0,
       y(2): 0,
       y(3): 0,
       y(4): 0,
       y(5): 0,
       y(6): 0,
       y(7): 0,
       y(8): 0,
       y(9): 0,
   }
)

但运行后得到的结果居然只是一个常数，完全没有依赖n的表达式，这到底是怎么回事？

原因分析 #

咱们先拆解一下问题，核心在于rsolve的能力边界：它本身更擅长处理低阶、形式标准的线性递推，比如常见的二阶、三阶递推。而你这个是10阶的递推，还带有跨10步的延迟项y(n-10)，这种复杂的结构已经超出了rsolve符号求解的舒适区。

再仔细看递推式的结构：非齐次项是常数0.04，rsolve在无法推导完整通解的情况下，可能直接返回了递推的**“伪稳态解”**——也就是假设当n趋向无穷大时，所有项都收敛到某个常数的情况。不过有意思的是，咱们手动代入验证的话，会发现这个递推其实不存在有限的稳态：假设稳态时$a_n = C$，代入递推式会得到$C = 0.996C + 0.004C + 0.04$，也就是$C = C + 0.04$，这显然矛盾，说明这个递推的解其实是随着n增长的，只是rsolve没能正确推导出来。

另外，你用sp.Rational转换小数的操作虽然没问题，但0.996转成分数是249/250，0.004是1/250，这种分数形式的递推特征方程是10次多项式，根的符号表达式极其复杂，Sympy根本无法直接求解并组合成通解，这也是它返回常数的原因之一。

解决办法 #

既然符号求解走不通，咱们可以换个思路：

• 改用数值计算：这是最直接的办法，用循环迭代计算每一项的数值，代码示例如下：

# 初始化前10项（n=0到9都是0）
a = [0.0 for _ in range(10)]
# 计算到n=100的结果，可按需调整范围
for i in range(10, 101):
    next_val = 0.996 * a[i-1] + 0.004 * a[i-10] + 0.04
    a.append(next_val)
# 打印n=10到n=15的结果看看
for idx, val in enumerate(a[10:16], start=10):
    print(f"a_{idx} = {val:.4f}")

• 尝试手动推导近似解：如果一定要符号形式，可以分析特征方程的根，比如其中一个根接近1，其他根的模小于1，这样通解可以近似为一个线性增长项加衰减项，但这种推导需要一定的数学功底，而且结果也会很复杂。
• 更新Sympy版本：旧版本的rsolve可能存在高阶递推处理的bug，更新到最新版说不定能解决问题，但别抱太大期望，毕竟10阶递推确实太复杂了。

备注：内容来源于stack exchange，提问作者Interpolated