I was trying to solve a non-linear system that consists of two equations with 2 unknowns, I tried using Sympy module to complete the task. After trying to run the code below (on Atom and Google Collab notebook), nothing is happening and the process is never terminated (waited for over that 15min). What could be the problem? is my system to complicated? I tried to use fsolve from scipy and i was able to get one of the roots in no time
This is the code below
import sympy as sym
sym.init_printing()
x,y = sym.symbols('x,y')
dmin = 70 #minimum of the actuator
dmax = 140 #maximum of the actuator
l3 = 10
l4 = 10
Amin = 0.523599 #min angle in radians
Amax = 2.0944 #max angle in radians
cos1 = sym.cos(Amin - sym.atan(l4/x) - sym.atan(l3/y))
cos2 = sym.cos(Amax - sym.atan(l4/x) - sym.atan(l3/y))
a2 = x**2 + l4**2
b2 = y**2 + l3**2
f=sym.Eq(a2 + b2 - 2*sym.sqrt(a2)*sym.sqrt(b2)*cos1, dmin**2)
g=sym.Eq(a2 + b2 - 2*sym.sqrt(a2)*sym.sqrt(b2)*cos2, dmax**2)
print(sym.solve([f,g],(x,y)))
Your equations are like:
In [160]: f
Out[160]:
__________ __________
2 2 ╱ 2 ╱ 2 ⎛ ⎛10⎞ ⎛10⎞ ⎞
x + y - 2⋅╲╱ x + 100 ⋅╲╱ y + 100 ⋅cos⎜atan⎜──⎟ + atan⎜──⎟ - 0.523599⎟ + 200 = 4900
⎝ ⎝x ⎠ ⎝y ⎠ ⎠
In [161]: g
Out[161]:
__________ __________
2 2 ╱ 2 ╱ 2 ⎛ ⎛10⎞ ⎛10⎞ ⎞
x + y - 2⋅╲╱ x + 100 ⋅╲╱ y + 100 ⋅cos⎜atan⎜──⎟ + atan⎜──⎟ - 2.0944⎟ + 200 = 19600
⎝ ⎝x ⎠ ⎝y ⎠ ⎠
Usually a mix of algebraic and transcendental functions like this would not have any anlytic solution but in this case the cos can cancel with the atan to make purely algebraic equations:
In [162]: expand_trig(f)
Out[162]:
__________ __________
2 2 ╱ 2 ╱ 2 ⎛ 0.866025291583566 5.00000194337561 5.00000194337561 86.6025291583566 ⎞
x + y - 2⋅╲╱ x + 100 ⋅╲╱ y + 100 ⋅⎜───────────────────────────── + ─────────────────────────────── + ─────────────────────────────── - ─────────────────────────────────⎟ + 200 = 4900
⎜ _________ _________ _________ _________ _________ _________ _________ _________⎟
⎜ ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ⎟
⎜ ╱ 1 + ─── ⋅ ╱ 1 + ─── y⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── x⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── x⋅y⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── ⎟
⎜ ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ⎟
⎝╲╱ x ╲╱ y ╲╱ x ╲╱ y ╲╱ x ╲╱ y ╲╱ x ╲╱ y ⎠
In [163]: expand_trig(g)
Out[163]:
__________ __________
2 2 ╱ 2 ╱ 2 ⎛ 0.500004241445914 8.6602295497065 8.6602295497065 50.0004241445914 ⎞
x + y - 2⋅╲╱ x + 100 ⋅╲╱ y + 100 ⋅⎜- ───────────────────────────── + ─────────────────────────────── + ─────────────────────────────── + ─────────────────────────────────⎟ + 200 = 19600
⎜ _________ _________ _________ _________ _________ _________ _________ _________⎟
⎜ ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ╱ 100 ⎟
⎜ ╱ 1 + ─── ⋅ ╱ 1 + ─── y⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── x⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── x⋅y⋅ ╱ 1 + ─── ⋅ ╱ 1 + ─── ⎟
⎜ ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ╱ 2 ⎟
⎝ ╲╱ x ╲╱ y ╲╱ x ╲╱ y ╲╱ x ╲╱ y ╲╱ x ╲╱ y ⎠
Using this we can massage the system into a pair of polynomials and then solve it with Groebner bases:
In [181]: from sympy.solvers.solvers import unrad
In [174]: eq1 = trigsimp(unrad(expand_trig(nsimplify(f)).rewrite(Add))[0]/(x*y)).factor().args[-1]
In [175]: eq2 = trigsimp(unrad(expand_trig(nsimplify(g)).rewrite(Add))[0]/(x*y)).factor().args[-1]
In [176]: gb = groebner([nsimplify(eq1.n()), nsimplify(eq2.n())], [x, y])
In [177]: for p in gb: print(p.n(2))
x - 2.9e-27*y**15 - 4.4e-25*y**14 + 2.8e-22*y**13 + 4.2e-20*y**12 - 1.2e-17*y**11 - 1.6e-15*y**10 + 2.7e-13*y**9 + 3.2e-11*y**8 - 3.6e-9*y**7 - 3.6e-7*y**6 + 2.7e-5*y**5 + 0.0021*y**4 - 0.099*y**3 - 6.2*y**2 + 1.3e+2*y + 6.5e+3
y**16 + 1.4e+2*y**15 - 1.0e+5*y**14 - 1.3e+7*y**13 + 4.4e+9*y**12 + 4.7e+11*y**11 - 1.1e+14*y**10 - 9.0e+15*y**9 + 1.5e+18*y**8 + 9.3e+19*y**7 - 1.3e+22*y**6 - 5.1e+23*y**5 + 6.3e+25*y**4 + 1.3e+27*y**3 - 1.5e+29*y**2 - 1.2e+30*y + 1.3e+32
In [178]: ry = [r.n() for r in real_roots(gb[-1])]
In [179]: ry
Out[179]:
[-112.210054906433, -94.6849042632114, -53.6834203472101, -49.4635109350851, 48.579017974275, 50.194835
2176397, 82.5120279645666, 103.812504924146, 118.173390238157, 142.710285016351]
In [180]: for yval in ry:
...: print('y =', yval, ', x =', solve(gb[0].subs(y, yval), x)[0])
...:
y = -112.210054906433 , x = 48.5790179742034
y = -94.6849042632114 , x = -53.6834203471953
y = -53.6834203472101 , x = -94.6849042632093
y = -49.4635109350851 , x = 103.812504924146
y = 48.5790179742750 , x = -112.210054906433
y = 50.1948352176397 , x = 118.173390238151
y = 82.5120279645666 , x = 142.710285016423
y = 103.812504924146 , x = -49.4635109348164
y = 118.173390238157 , x = 50.1948352170293
y = 142.710285016351 , x = 82.5120279639959
Rewriting as polynomial introduces spurious solutions so only some of those are actually solutions of f and g.
A simpler and faster way to compute these solutions if you have an initial guess for a particular solution that you want is just to use nsolve to solve the system numerically:
In [194]: nsolve([f, g], [x, y], [118, 50])
Out[194]:
⎡118.173390238157⎤
⎢ ⎥
⎣50.1948352176396⎦