Symbolic manipulation of equation with substitutions given relations between variables in Sympy

Question

I'm trying to use Python SymPy to change equations 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)

and 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))

into equations 2

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

and 4

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

respectively using the substitutions

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)

I attached an image of the equations and substitutions.

Equations and substitutions

I tried the below code for energy and momentum, which solves for F². (I ultimately want the reciprocal 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

But the output for equation 1 (energy) does not match equation 3. There are still some lowercase letters in the output for equation 3 (momentum), which is not what I want. I attached images of manipulating by hand the equations 1 (energy) and 3 (momentum) into equations 2 and 4, which is what I want Python to do for me instead of calculating manually.

Manipulation by hand for energy (equations 1 and 3)

Manipulation by hand for momentum (equations 2 and 4)

Could someone please enlighten me? Thank you very much!

Answer 1

As far as I know, the only way to achieve the results you are looking for with SymPy is by following step by step what you would do by hand. There is no magic function where you can provide your original equation with a list of substitution and expect the correct form to come out. The search space is just too big.

However, I assume you would like to see if SymPy is actually any good at expression manipulation. The short answer is: it depends on the task at hand. I'm going to provide you with only the first part of the manipulation of your energy equation. The process is long and tedious and I run out of time, so I only achieved a partial "result". What is missing is for you to figure it out, gaining experience by practice :)

I'm going to use Algebra With SymPy because it provides a nice Equation class which can be used with mathematical operations. Depending on the task at hand, it could make life a little bit easier...

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

Let's start by diving by rho1:

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

Expanding the terms in order to get rho2/rho1 is easy. The difficult part would be to re-collect the terms in a logical way. Instead, I got Req2 with a simple hand manipulation, which saves some typing and troubles:

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

Note the second term on the left hand side. We need to select it, expand it and manipulate it further:

# `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

Now, let's do the rest:

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

Note that sympy didn't perform the trivial simplification (R - 1) / (1 - R) = -1. We can apply like this:

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

I've run out of time, the rest is up to you. As you can see, the process is very much "by hand". So, on top of the difficulties of doing it by hand, you are also increasing the difficulty even more by bumping into issues that are trivial to solve with hand-manipulation, but not so trivial with sympy....