python - 我无法将线性方程符号化地求解为变量，即使手动完成它非常简单

Question

我目前非常沮丧，因为我有一个非常简单的方程式，其中的符号无法求解为变量。即使手动完成，它10^{-17}在替换时也会出现错误（数字上为 的顺序）。

我无法发布 mwe 示例，因此我将发布问题的最简化形式。

import sympy as sy
from sympy import sqrt, Derivative
r          = sy.Symbol('r')
uf         = sy.Function('uf')
Expression =-4.5*r*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) - (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62))*(r**2 - 2*r + 0.81)*uf(r)/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*sqrt(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r - (r**2 + 3.03701355989026e-33)*uf(r)**2 + 1.0*(r**2 + 0.81)**2 - 0.6561)) - 1.0*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(r**2 - 2*r + 0.81)*(1.0*r*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(r**2 - 2*r + 0.81)*uf(r)/((r**2 + 3.03701355989026e-33)**(3/2)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)) + 0.5*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62))*(r**2 - 2*r + 0.81)*uf(r)/((r**2 + 3.03701355989026e-33)**(3/2)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)**2) + 0.5*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(-2*r*(r**2 - 2*r + 0.81)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*uf(r) + 4*(r - 1)*(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*uf(r) + (r**2 + 3.03701355989026e-33)*(2*(1 - r)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 - 2*r + 0.81)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62))*uf(r) + 2*(r**2 + 3.03701355989026e-33)*(r**2 - 2*r + 0.81)*(-4*r*(r**2 + 0.81) + 1.62*r - 1.62)*uf(r) + 2*(r**2 + 3.03701355989026e-33)*(r**2 - 2*r + 0.81)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*Derivative(uf(r), r) + 2*(r*(r**2 - 2*r + 0.81) + (1 - r)*(r**2 + 3.03701355989026e-33))*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*uf(r))/((r**2 + 3.03701355989026e-33)**(3/2)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)**2))*uf(r)/(sqrt(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)) + 0.5*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(r**2 - 2*r + 0.81)*(-1.8*r*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 + 3.03701355989026e-33)**(7/2)*(2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) - (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62))*(r**2 - 2*r + 0.81)*sqrt(1/(r**2 - 2*r + 0.81))*uf(r) - 2.0*r*(-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))**2*(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)**2 + 1.0*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 + 3.03701355989026e-33)**(7/2)*(1.8*r**2 - 5.46662440780247e-33)*(r**2 - 2*r + 0.81)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*sqrt(1/(r**2 - 2*r + 0.81))*uf(r) + 1.0*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 + 3.03701355989026e-33)**(7/2)*(1.8*r*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62)) + (1.8*r**2 - 5.46662440780247e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 - 2*r + 0.81)*sqrt(1/(r**2 - 2*r + 0.81))*uf(r) + 2.0*sqrt((-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)/(r**2 - 2*r + 0.81))*(r**2 + 3.03701355989026e-33)**4*(-0.9*r*(1.8*r*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62)) + (1.8*r**2 - 5.46662440780247e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)) - (0.9*r*(1.8*r**2 - 5.46662440780247e-33) - (r**2 - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 - 2*r + 0.81)**2*sqrt(1/(r**2 - 2*r + 0.81)))/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))**2*(r**2 + 3.03701355989026e-33)**5*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)**2*sqrt(1/(r**2 - 2*r + 0.81))) + 2.5*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(1.8*r**2 - 5.46662440780247e-33)*(r**2 - 2*r + 0.81)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*uf(r)/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*sqrt(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r - (r**2 + 3.03701355989026e-33)*uf(r)**2 + 1.0*(r**2 + 0.81)**2 - 0.6561)) + 2.5*sqrt((r**2 - 2*r + 0.81)/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(1.8*r*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62)) + (1.8*r**2 - 5.46662440780247e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 - 2*r + 0.81)*uf(r)/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*sqrt(r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r - (r**2 + 3.03701355989026e-33)*uf(r)**2 + 1.0*(r**2 + 0.81)**2 - 0.6561)) - 2.5*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62))*(r**2 - 2*r + 0.81)*uf(r)**2/((r**2 + 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)*(-0.81*r**2 + 1.62*r - (r**2 + 3.03701355989026e-33)*uf(r)**2 + 1.0*(r**2 + 0.81)**2 - 0.6561)) + (-0.9*r*(r**2 + 3.03701355989026e-33)*(1.8*r*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62)) + (1.8*r**2 - 5.46662440780247e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 - 2*r + 0.81)/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))**2*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)) - (r**2 + 3.03701355989026e-33)*(0.9*r*(1.8*r**2 - 5.46662440780247e-33) - (r**2 - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(r**2 - 2*r + 0.81)/(-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))**2)*(r**2 - 2*r + 0.81)/(r**2 + 3.03701355989026e-33) + (-0.9*r*(1.8*r*(-2*r*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + (r**2 + 3.03701355989026e-33)*(4*r*(r**2 + 0.81) - 1.62*r + 1.62)) + (1.8*r**2 - 5.46662440780247e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561)) - (0.9*r*(1.8*r**2 - 5.46662440780247e-33) - (r**2 - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))*(1.0*(r**2 + 3.03701355989026e-33)*uf(r)**2/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + 1)*(r**2 - 2*r + 0.81)**2/((-3.24*r**2 + (-r**2 + 2*r - 3.03701355989026e-33)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))**2*(-1.0*(r**2 + 3.03701355989026e-33)*uf(r)**2/(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561) + 1.0)*(-0.81*r**2 + 1.62*r + (r**2 + 0.81)**2 - 0.6561))
sy.solve(Expression,sy.diff(uf(r),r)) 

您可以验证方程很容易求解，sy.diff(uf(r),r)因为 sy.diff(Expression,sy.diff(uf(r),r),1)它是非零并且sy.diff(Expression,sy.diff(uf(r),r),2)是0。

上面的求解器在使用 jupyter 时会破坏我的内核，我必须手动求解solution = (Expression-sy.diff(Expression,sy.diff(uf(r),r),1)*sy.diff(uf(r),r))/(-sy.diff(Expression,sy.diff(uf(r),r),1))——仍然不完全准确

标签： pythonmathsympy

Answer 1

即使这是一个笨拙的表达式，也可以使用 SymPy 工具解决这个问题。目标是减少必须跟踪和操作的符号数量：

0) 将导数替换为 z 以便于跟踪

>>> from sympy import diff
>>> from sympy.abc import z
>>> y = Expression.subs(diff(uf(r), r), z)

1) 提取常用子表达式

>>> from sympy import cse
>>> reps, e = cse(y)

2) e 是一个列表，得到表达式

>>> e = e[0]

3) 展开 e

>>> e = e.expand(); type(e)
Add

4）隔离z

>>> i, d = e.as_independent(z); type(d)
Mul

5) 求解 z

>>> sol = i/d.coeff(z)

6）恢复提取的模式（以相反的顺序）

>>> sol = sol.subs(reps[::-1]); sol.count_ops()
1584

把它们放在一起，我们有这个后备：

def unwieldy_linear_solve(eq, x):
    from sympy import cse
    if not x.is_Symbol:
        d = Dummy()
        return unwieldy_linear_solve(eq.subs(x, d), d)
    r, e = cse(eq)
    i, d = e[0].expand().as_independent(x)
    assert d == x or d.is_Mul and x in d.args
    c, z = d.as_independent(x)
    assert z == x
    sol = i/c
    for o,n in r[::-1]:
        sol = sol.xreplace({o:n})
    return sol

>>> count_ops(unwieldy_linear_solve(Expression, diff(uf(r),r)))
1584