sympy: solving for functions results in terms of other functions

Question

I have this system of equations in sympy:

m1, m2 = symbols('m1, m2', real=True, positive=True)
t0, t = symbols('t0 t', real=True, positive=True)
p1, p2, v1, v2, a1, a2, f1, f2 = symbols('p1 p2 v1 v2 a1 a2 f1 f2', cls=Function)

G = 6.67384e-11

eq1 = Eq(f1(t), G * m1 * m2 * (p2(t) - p1(t)) / abs(p2(t) - p1(t)))
eq2 = Eq(a1(t), a1(t0) + (1/m1) * Integral(f1(t), (t, t0, t)))
eq3 = Eq(v1(t), v1(t0) + Integral(a1(t), (t, t0, t)))
eq4 = Eq(p1(t), p1(t0) + Integral(v1(t), (t, t0, t)))

eq5 = Eq(f2(t), G * m1 * m2 * (p1(t) - p2(t)) / abs(p2(t) - p1(t)))
eq6 = Eq(a2(t), a2(t0) + (1/m2) * Integral(f2(t), (t, t0, t)))
eq7 = Eq(v2(t), v2(t0) + Integral(a2(t), (t, t0, t)))
eq8 = Eq(p2(t), p2(t0) + Integral(v2(t), (t, t0, t)))

I want to have sympy solve these equations for p1(t) and p2(t) in terms of t, but if I use this:

results = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8],
                [p1(t), p2(t)])

I simply get p1(t) and p2(t2) in terms of v1(t) and v2(t):

p1(t0) + Integral(v1(t), (t, t0, t))
p2(t0) + Integral(v2(t), (t, t0, t))

If I use this:

results = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8],
                [p1(t), p2(t), v1(t), v2(t), a1(t), a2(t), f1(t), f2(t)])

I get:

p1(t0) + Integral(v1(t0) + Integral(a1(t), (t, t0, t)), (t, t0, t))
p2(t0) + Integral(v2(t0) + Integral(a2(t), (t, t0, t)), (t, t0, t))

(I've omitted the results for the other functions because they're very long.)

How can I make sympy to calculate the results for p1(t) and p2(t) in pure terms of t?

Answer 1

Try replacing (p2(t) - p1(t)) / abs(p2(t) - p1(t))) with Symbol('s') and (p1(t) - p2(t)) / abs(p2(t) - p1(t))) with -Symbol('s'). Then re-solve for the variables of interest, backsubstituting as necessary:

>>> from sympy.utilities.misc import filldedent
>>> def backsub(sol):
...     for i in range(len(sol)):
...         was = sol.copy()
...         for k in sol.keys():
...             sol[k] = sol[k].subs(sol)
...         if was == sol:
...             break
...
>>> v
(p1, p2, v1, v2, a1, a2, f1, f2)
>>> v=[i(t) for i in _]
>>> sol=solve(eqs,v);backsub(sol)
>>> print filldedent(sol[p1(t)].doit().simplify())

1.11230666666667e-11*m2*s*t**3 - 3.33692e-11*m2*s*t**2*t0 +
3.33692e-11*m2*s*t*t0**2 - 1.11230666666667e-11*m2*s*t0**3 +
0.5*t**2*a1(t0) - 1.0*t*t0*a1(t0) + 1.0*t*v1(t0) + 0.5*t0**2*a1(t0) -
1.0*t0*v1(t0) + 1.0*p1(t0)
>>> print filldedent(sol[p2(t)].doit().simplify())

-1.11230666666667e-11*m1*s*t**3 + 3.33692e-11*m1*s*t**2*t0 -
3.33692e-11*m1*s*t*t0**2 + 1.11230666666667e-11*m1*s*t0**3 +
0.5*t**2*a2(t0) - 1.0*t*t0*a2(t0) + 1.0*t*v2(t0) + 0.5*t0**2*a2(t0) -
1.0*t0*v2(t0) + 1.0*p2(t0)