在 Sympy 中给定变量之间的关系，通过替换对方程进行符号操作

Question

我正在尝试使用 Python SymPy 来更改方程 1

#Equation 1
energy = sp.Eq(sp.together((rho2*c2prime**2-rho1*c1prime**2)/ 2)
                + g*(rho2-rho1)*h1prime,
               sp.together((rho2*c**2-rho1*c**2)/ 2)
                + g*(rho2-rho1)*h1)

和3

#Equation 3
momentum = sp.Eq((rho1*h1+rho2*h2)*c**2,
                 sp.Rational(1,2)*rho1*h1prime*(c**2+c1prime**2)
                  + sp.Rational(1,2)*rho2*h2prime*(c**2+c2prime**2)
                  - sp.Rational(1,2)*rho1*g*(h1-h1prime)**2
                  - sp.Rational(1,2)*rho2*g*(h2-h2prime)**2
                  - rho2*g*(h1-h1prime)*(h2-h2prime))

代入方程2

#Equation 2
energy7 = sp.Eq(1/F_sq,(H1+H1prime)/(2*H1prime**2)
                      +(((2-(H1+H1prime))*R)/(2*(1-H1prime)**2)))

和4

#Equation 4
momentum6 = sp.Eq(1/F_sq,1/H1prime+R/(1-H1prime))

分别使用替换

R=rho2/rho1
H1=h1/(h1+h2)
H2=h2/(h1+h2)
H1prime=h1prime/(h1+h2)
H2prime=h2prime/(h1+h2)
c1prime=c*H1/H1prime
c2prime=c*H2/H2prime
F_sq=c**2/(g*(h1+h2)*(1-rho2/rho1))
H2=1-H1
H2prime=1-H1prime
H2-H2prime=H1prime-H1
H1-H1prime=H2prime-H2
H2prime+H2=2-(H1prime+H1)

我附上了方程式和替换的图像。

我尝试了下面的能量和动量代码，它求解了 F²。 （我最终想要倒数 1/F²。）

energy2 = energy.subs({c1prime:c*H1/H1prime, c2prime:c*H2/H2prime})
energy3 = energy2.subs({h1:H1*(h1+h2), h1prime:H1prime*(h1+h2)})
energy4 = energy3.subs({rho2:R*rho1})
energy5 = energy4.subs({h1+h2:c**2/(g*F_sq*(1-R))})
energy6 = energy5.subs({H2prime:1-H1prime})
energy7 = sp.solve(energy6, F_sq)[0].simplify() #Not Equation 2

momentum2 = momentum.subs({c1prime:c*H1/H1prime, c2prime:c*H2/H2prime})
momentum3 = momentum2.subs({h1:H1*(h1+h2), h1prime:H1prime*(h1+h2)})
momentum4 = momentum3.subs({rho2:R*rho1})
momentum5 = momentum4.subs({h1+h2:c**2/(g*F_sq*(1-R))})
momentum6 = sp.solve(momentum5, F_sq)[0].simplify() #Not Equation 4

但是方程1（能量）的输出与方程3不匹配。方程3（动量）的输出中仍然有一些小写字母，这不是我想要的。我将手动操作方程 1（能量）和 3（动量）的图像附加到方程 2 和 4 中，这就是我希望 Python 为我做的事情，而不是手动计算。

有人可以启发我吗？非常感谢！

Answer 1

据我所知，使用 SymPy 实现您想要的结果的唯一方法是逐步遵循您手动执行的操作。没有什么神奇的功能可以让您为原始方程提供替换列表并期望得到正确的形式。搜索空间太大了。

但是，我假设您想看看 SymPy 是否真的擅长表达操作。简短的回答是：这取决于手头的任务。我将只向您提供能量方程操作的第一部分。这个过程漫长而乏味，而且我的时间也不够了，所以我只取得了部分“结果”。缺少的是你要弄清楚，通过实践获得经验:)

我将使用 Algebra With SymPy，因为它提供了一个很好的

Equation

类，可以与数学运算一起使用。根据手头的任务，它可以让生活变得更轻松......

from sympy import *
from algebra_with_sympy import Eqn, algwsym_config, solve as esolve
algwsym_config.output.solve_to_list = True
algwsym_config.output.label = False

c, g, R, F = symbols("c, g, R, F")
c1, c2, c1prime, c2prime = symbols(r"c1, c2, {c_1}', {c_2}'")
h1, h2, h1prime, h2prime = symbols("h1, h2, {h_1}', {h_2}'")
H1, H2, H1prime, H2prime = symbols("H1, H2, {H_1}', {H_2}'")
rho1, rho2, r1prime, r2prime = symbols("rho1, rho2, {rho_1}', {rho_2}'")

# use the linearity property in order to apply a function representing
# some operation to each term of an additive expression
def apply_term_by_term(expr, function):
    if isinstance(expr, (Eq, Eqn)):
        return expr.func(*[apply_term_by_term(side, function) for side in expr.args])
    if isinstance(expr, Add):
        return expr.func(*[function(t) for t in expr.args])
    return expr

energy = Eqn(
    together((rho2*c2prime**2-rho1*c1prime**2)/ 2) + g*(rho2-rho1)*h1prime,
    together((rho2*c**2-rho1*c**2)/ 2) + g*(rho2-rho1)*h1)
energy

让我们从

rho1

潜水开始：

energy = apply_term_by_term(energy, lambda t: t / rho1)
energy

扩展术语以获得

rho2/rho1

很容易。困难的部分是以逻辑方式重新收集术语。相反，我通过简单的手动操作得到了

Req2

，这节省了一些打字和麻烦：

Req1 = Eqn(R, rho2 / rho1)
Req2 = Eqn((rho2 - rho1) / rho1, R - 1)
energy = energy.subs(Req2)
energy

注意左侧的第二项。我们需要选择它，展开它并进一步操作它：

# `find` returns a set of results. When a set is converted
# to a list, the order of elements might change from run
# to run. It is not reproducible. One way to get reproducible
# results is to somehow sort the list of results.
# I'm going to find a set of multiplications containing `c1prime`,
# sorting them according to the number of operations.
# The last element of the list (the bigger term) will be the
# term I'm looking for.
second_term_on_lhs = sorted(
    energy.lhs.find(lambda t: isinstance(t, Mul) and t.has(c1prime)),
    key=count_ops
)[1]
new_second_term_on_lhs = second_term_on_lhs.expand().subs(*Req1.swap.args).together()
new_second_term_on_lhs

energy = energy.subs(second_term_on_lhs, new_second_term_on_lhs)
energy

现在，让我们做剩下的事情：

energy = apply_term_by_term(energy, lambda t: t / (g*(h1 + h2)))
energy

H1eq = Eqn(H1, h1 / (h1+h2))
H1peq = Eqn(H1prime, h1prime / (h1+h2))
energy = energy.subs(H1eq.swap, H1peq.swap)
energy

energy = apply_term_by_term(energy, lambda t: t / (1 - R))
energy

请注意，sympy 没有执行简单的简化

(R - 1) / (1 - R) = -1

。我们可以这样申请：

energy = apply_term_by_term(asd, lambda t: t.simplify() if t.has(R - 1) else t)
energy

c1prime_eq = Eqn(c1prime, c*H1/H1prime)
c2prime_eq = Eqn(c2prime, c*H2/H2prime)
energy = energy.subs(c1prime_eq, c2prime_eq).collect(c)
energy

F_eq = Eqn(F**2, c**2 / (g * (h1 + h2) * (1 - R)))
energy = energy.subs(F_eq.swap)
energy

second_term_on_rhs = sorted(
    energy.rhs.find(lambda t: isinstance(t, Mul) and t.has(c)),
    key=count_ops
)[0]
new_second_term_on_rhs = (second_term_on_rhs / (1 - R)).subs(*F_eq.swap.args) * (1 - R)
new_second_term_on_rhs

energy = energy.subs(second_term_on_rhs, new_second_term_on_rhs)
energy

energy = energy - new_second_term_on_rhs + H1prime
energy

我的时间不多了，剩下的就看你的了。正如您所看到的，该过程很大程度上是“手动”的。因此，除了手动执行此操作的困难之外，您还会遇到通过手动操作解决的微不足道的问题，但通过 sympy 解决却不是那么微不足道的问题，从而进一步增加了难度......