python - Sympy 与 Cython 在尝试查找根时生成奇异矩阵

Question

我正在尝试使用 Sympy 求解一组非线性方程。这是我的代码，更改了一些数字，并且只有一个输入值(真正的代码运行了 170,000 行数据):

import sympy as sp

K0 = 2.00*10**-4
x2, y2 = sp.symbols('x2, y2')
x0 = 500
y0 = 500
R1 = ((x2 - x0)**2 + (y2 - y0)**2)**0.5
R2 = K0 * R1
H2 = sp.atan(R2 * (x2 - x0)/R1)
V2 = sp.atan(R2 * (y2 - y0)/R1)
x, y = sp.symbols('x, y')
x0 = 1.0
y0 = 2.0
x = R1 * H2
y = R1 * V2
dat = sp.nsolve([x - x0, y - y0], [x2, y2], [512, 512]) # This line is the problem
print "dat = %f, %f" % (dat[0], dat[1])

纯粹使用 Python，代码运行良好并产生良好的输出 - 但速度非常慢。为了加快这个过程，我使用 Cython 编译了一个具有完全相同代码的 .pyx 文件(在定义中，def test())，但在运行时我遇到了:

  File "test.py", line 3, in <module>
    demo.test()
  File "demo.pyx", line 17, in demo.test
    dat = sp.nsolve([x - x0, y - y0], [x2, y2], [512, 512])
  File "C:\...\site-packages\sympy\utilities\decorator.py", line 91, in func_wrapper
    return func(*args, **kwargs)
  File "C:\...\site-packages\sympy\solvers\solvers.py", line 2847, in nsolve
    x = findroot(f, x0, J=J, **kwargs)
  File "C:\...\site-packages\mpmath\calculus\optimization.py", line 960, in findroot
    for x, error in iterations:
  File "C:\...\site-packages\mpmath\calculus\optimization.py", line 658, in __iter__
    s = self.ctx.lu_solve(Jx, fxn)
  File "C:\...\site-packages\mpmath\matrices\linalg.py", line 227, in lu_solve
    A, p = ctx.LU_decomp(A)
  File "C:\...\site-packages\mpmath\matrices\linalg.py", line 137, in LU_decomp
    raise ZeroDivisionError('matrix is numerically singular')
ZeroDivisionError: matrix is numerically singular

我已将问题范围缩小到 x - x0 和 y - y0 部分。由于某种原因，编译后的代码无法处理不等于 0 时求根的问题。 nsolve 是否可以简单地使用 Cython 转换为 C 语言？我是否遗漏了一些与 Cython 有关的内容？

Answer 1

您可以将sympy.lambdify与例如一起使用SciPy 的求解器。如果这还不够快，您可以使用symengine.Lambdify。

获得正确的函数签名并导出雅可比行列式需要您克服重重困难。如果你想使用我写过的库 pyneqsys :

>>> from pyneqsys.symbolic import SymbolicSys
>>> neqsys = SymbolicSys([x2, y2], [x - x0, y - y0])
>>> neqsys.solve([512, 512])
Out[4]: 
(array([ 547.28609349,  594.58064617]),
     fjac: array([[ 0.91320338,  0.4075041 ],
       [-0.4075041 ,  0.91320338]])
     fun: array([ -1.37667655e-13,   1.52011737e-12])
 message: 'The solution converged.'
    nfev: 17
    njev: 2
     qtf: array([  1.55620322e-10,   4.63225371e-10])
       r: array([ 0.02751454,  0.023682  ,  0.03261983])
  status: 1
 success: True
       x: array([ 547.28609349,  594.58064617]))

如果这 170 000 个解决方案涉及逐渐改变参数，pyneqsys 可以利用这一点(通过在解决方案之间将解决方案作为猜测传播)。它还可以通过设置环境变量SYM_BACKEND=sympysymengine

自动使用 symengine.Lambdify