为什么 SymPy 的 linsolve 函数无法求解此问题中倒数第二个变量？答案

【问题标题】：Why is SymPy's linsolve function not able to solve for the second to the last variable in this problem?为什么 SymPy 的 linsolve 函数无法求解此问题中倒数第二个变量？
【发布时间】：2021-08-02 06:11:03
【问题描述】：

(TLDR：当您将方程作为纯 Python 列表传递时，SymPy 的 linsolve 函数无法求解线性方程组，该方程组是通过将有限差分法应用于 ODE BVP 生成的，但能够这样做说SymPy 的 Matrix 函数中的方程式列表。这可能是一个需要修复的错误，特别是考虑到在 SymPy 文档中，the example 他们给你的让他们将列表作为参数传递给 linsolve。）

我有一个边值问题常微分方程，我打算使用有限差分法求解。我的 ODE，在 SymPy 表示中，是x*y(x).diff(x,2) + y(x).diff(x) + 500 = 0，有y(1)=600 和y(3.5)=25。整个代码如下：

import sympy as sp
from numpy import *
from matplotlib.pyplot import *
y = sp.Function('y')
ti = time.time()
steps = 10
ys = [ y(i) for i in range(steps+1) ]
ys[0], ys[-1] = 600, 25
xs = linspace(1, 3.5, steps+1)
dx = xs[1] - xs[0]
eqns = [ xs[i]*(ys[i-1] - 2*ys[i] + ys[i+1])/dx**2 + 
        ( ys[i+1] - ys[i-1] )/2/dx + 500
        for i in range(1, steps) ]
ys[1:-1] = sp.linsolve(eqns, ys[1:-1]).args[0]
scatter(xs, ys)
tf = time.time()
print(f'Time elapsed: {tf-ti}')

十个步骤，这工作得很好。但是，如果我立即超过 11 步，SymPy 将不再能够求解线性方程组。尝试绘制结果会引发TypeError: can't convert expression to float。检查 y 值 ys 的列表会发现其中一个变量，特别是倒数第二个变量，并没有被 sympy 的 linsolve 解决，而是根据这个未解决的变量来解决其他变量的解决方案。例如，如果我有 50 个步骤，倒数第二个变量是 y(49)，这在其他未知数的解决方案中具有特色，并且未解决，例如2.02061855670103*y(49) - 26.1340206185567。相比之下，在我解决的另一个 BVP ODE 中，y(x).diff(x,2) + y(x).diff(x) + y(x) - 1 与 y(0)=1.5 和 y(3)=2.5，无论我想要 10、50 还是 200 步都没有任何问题。它很好地解决了所有变量，但这似乎是一个特殊的例外，因为我在许多其他 ODE 中遇到了上述问题。

SymPy 的不一致在这里非常令人沮丧。唯一的安慰是，在我遇到这个问题之前，我实际上已经用我想要的不同数量的步骤解决了几次。我将eqns 变量包含在sympy 的Matrix 函数中，就像ys[1:-1] = sp.linsolve(sp.Matrix(eqns), ys[1:-1]).args[0] 一样，只是因为它在终端中以这种方式显示得更好。但为了解决在脚本文件中，我认为将其包裹在 sp.Matrix 中是不必要的，我自然将其删除以简化事情。

【问题讨论】：

标签： python-3.x matrix sympy differential-equations finite-difference

【解决方案1】：

为 SO（或其他任何地方）格式化问题以提供完整的代码示例而不会丢失导入等是有礼貌的。此外，您应该将问题精简到最小情况并删除所有不必要的细节。考虑到这一点，证明问题的更好方法是实际弄清楚 linsolve 的参数是什么并直接呈现它们，例如：

from sympy import *

y = Function('y')

eqns = [
    -47.52*y(1) + 25.96*y(2) + 13436.0,
    25.96*y(1) - 56.32*y(2) + 30.36*y(3) + 500,
    30.36*y(2) - 65.12*y(3) + 34.76*y(4) + 500,
    34.76*y(3) - 73.92*y(4) + 39.16*y(5) + 500,
    39.16*y(4) - 82.72*y(5) + 43.56*y(6) + 500,
    43.56*y(5) - 91.52*y(6) + 47.96*y(7) + 500,
    47.96*y(6) - 100.32*y(7) + 52.36*y(8) + 500,
    52.36*y(7) - 109.12*y(8) + 56.76*y(9) + 500,
    56.76*y(8) - 117.92*y(9) + 61.16*y(10) + 500,
    61.16*y(9) - 126.72*y(10) + 2139.0
]

syms = [y(1), y(2), y(3), y(4), y(5), y(6), y(7), y(8), y(9), y(10)]

print(linsolve(eqns, syms))

您希望在这里得到每个未知数的简单数值解，但返回的结果（来自 SymPy 1.8）是：

FiniteSet((5.88050359812056*y(10) - 5.77315239260531, 10.7643116711359*y(10) - 528.13328974178, 14.9403214726998*y(10) - 991.258097567359, 9.85496358613721e+15*y(10) - 1.00932650309452e+18, 7.35110502818395*y(10) - 312.312287998229, 5.84605452313345*y(10) - 217.293922525318, 4.47908204606922*y(10) - 141.418192750506, 3.22698120573309*y(10) - 81.4678489766327, 2.07194244604317*y(10) - 34.9738391105298, 1.0*y(10)))

如果系统没有唯一解，linsolve 函数将返回一个包含一个或多个未知数的解。另请注意，有一些像9.85496358613721e+15 这样的大数字表明这里可能存在数字问题。

其实这是 SymPy 中的一个 bug，它已经在 master 分支上修复了： https://github.com/sympy/sympy/pull/21527

如果您从 git 安装 SymPy，那么您可以找到以下输出：

FiniteSet((596.496767861074, 574.326903264955, 538.901024315178, 493.575084012669, 440.573815245681, 381.447789421181, 317.320815173574, 249.033388036155, 177.23053946471, 102.418085492911))

还请注意，通常最好避免在 SymPy 中使用浮点数，因为它是一个专为精确符号计算而设计的库。使用 NumPy 或其他一些可以使用 BLAS/LAPACK 例程的固定精度浮点库，可以更有效地求解这样的浮点方程系统。要在此处使用带 sympy 的有理算术，您只需将 linspace 行更改为

xs = [sp.Rational(x) for x in linspace(1, 3.5, steps+1)]

它可以在 SymPy 1.8 上正常工作，实际上速度更快（至少如果你安装了 gmpy2）。

【讨论】：