sympy hangs when trying to solve for simultaneous equations with cubics and trig

Question

I'm trying to solve for this question: https://math.stackexchange.com/questions/4370530/equation-of-a-cubic-bezier-curve/4370920#4370920

Here's the full code

from sympy import *

xa=symbols('xa')
xb=symbols('xb')
xc=symbols('xc')
xd=symbols('xd')
xe=symbols('xe')
xf=symbols('xf')
xg=symbols('xg')
xh=symbols('xh')
xi=symbols('xi')
xm=symbols('xm')
xn=symbols('xn')
ya=symbols('ya')
yb=symbols('yb')
yc=symbols('yc')
yd=symbols('yd')
ye=symbols('ye')
yf=symbols('yf')
yg=symbols('yg')
yh=symbols('yh')
yi=symbols('yi')
ym=symbols('ym')
yn=symbols('yn')

p=symbols('p')
q=symbols('q')

alpha,beta=symbols("alpha,beta")

pxe = xa+p*(xb-xa)
pxf = xb+p*(xc-xb)
pxg = xc+p*(xd-xc)
pxh = xe+p*(xf-xe)
pxh = pxh.subs(xe, pxe).subs(xf, pxf)
pxi = xf+p*(xg-xf)
pxi = pxi.subs(xf, pxf).subs(xg, pxg)
pxm = xh+p*(xi-xh)
pxm = pxm.subs(xh, pxh).subs(xi, pxi)
pxm = simplify(pxm)
pym = pxm.subs(xa, ya).subs(xb,yb).subs(xc,yc).subs(xd,yd)

ab = symbols("ab")
cd = symbols("cd")
lxb = xa + ab*sin(alpha)
lxc = xd + cd*sin(beta)
lyb = ya + ab*cos(alpha)
lyc = yd + cd*cos(beta)

pxm = pxm.subs(xb,lxb).subs(xc,lxc)
pym = pym.subs(yb,lyb).subs(yc,lyc)

eqxm = Eq(xm, pxm)
cdx = solve(eqxm, cd)[0]
eqym = Eq(ym,pym)
cdy = solve(eqym, cd)[0]
eqcd = Eq(cdx, cdy)
abp = solve(eqcd, ab)[0]
abq = abp.subs(p,q).subs(xm, xn).subs(ym,yn)
print(abp)
print(abq)
eqpq = Eq(abp, abq)
qv = solve(eqpq, q)
print(qv)

I get

abp = -(2*p**3*xa*cos(beta) - 2*p**3*xd*cos(beta) - 2*p**3*ya*sin(beta) + 2*p**3*yd*sin(beta) - 3*p**2*xa*cos(beta) + 3*p**2*xd*cos(beta) + 3*p**2*ya*sin(beta) - 3*p**2*yd*sin(beta) + xa*cos(beta) - xm*cos(beta) - ya*sin(beta) + ym*sin(beta))/(3*p*(p - 1)**2*sin(alpha - beta))

abq = -(2*q**3*xa*cos(beta) - 2*q**3*xd*cos(beta) - 2*q**3*ya*sin(beta) + 2*q**3*yd*sin(beta) - 3*q**2*xa*cos(beta) + 3*q**2*xd*cos(beta) + 3*q**2*ya*sin(beta) - 3*q**2*yd*sin(beta) + xa*cos(beta) - xn*cos(beta) - ya*sin(beta) + yn*sin(beta))/(3*q*(q - 1)**2*sin(alpha - beta))

when I try to solve for q when abp = abq, python process hangs. Could anyone inform me the kinds of equations sympy has difficulty solving ?

Answer 1

The solution is unwieldy (it would probably be better to substitute in values for all but q and then use nsolve to solve the equation) but if you use simplify=False, check=False the solution will be obtained fairly quickly.

>>> gsol = solve(eqpq, q, check=False, simplify=False)
>>> reps = {xd: 10, beta: 11, ya: 12, alpha: 13,
...         p: 14, xn: 15, yn: 16, xa: 17, yd: 18, xm: 20, ym: 21}
>>> [i.n(3, subs=reps, chop=True) for i in gsol]  # slow
[0.416, 1.61, 13.9]
>>> [i.subs(reps).n(3,chop=True) for i in gsol]  # faster
[0.416, 1.61, 13.9]
>>> neq = eqpq.rewrite(Add).as_numer_denom()[0].subs(reps).n()
>>> [i.n(3) for i in real_roots(neq)]  # very fast
[0.416, 1.61, 13.9]

In the slowest case, more care is taken to make sure the answer is correct. In the last case, the numerical values are substituted into the expression and just the numerator is retained and the real roots of that are requested.