0

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

我有一个边值问题常微分方程，我打算使用有限差分法求解。在 SymPy 表示中，我的 ODE 是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是不必要的，我自然将其删除以简化事情。

4

1 回答 1

3

为 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 中的一个错误，它已经在主分支上修复： 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）。

于 2021-08-02T12:03:13.203 回答