Solving system of trig equations using sympy-CodePudding

I made a system of equations for an optimization of a pro-link motorcycle into matlab that I'm turning into python code beacuse i need to load it into another software. The matlab code is the following:

clc
clear 

Lmono= 320;
Lbielletta=145;
IPS= 16.03;
Pivot= -20;

R1x=(-46.72)-((Pivot)*sind(IPS));
R1y=((180.26)-((Pivot)*cosd(IPS)));
R_1x=(-43.52)-((Pivot)*sind(IPS));
R_1y=((-151.37)-((Pivot)*cosd(IPS)));
R4=60; 
R3=203.727;
Ip=100.1;
eta=36.79;


syms phi o 
eqns = [(Lbielletta)^2 == ((R3*cosd(phi)-R4*sind(o) R_1x)^2) ((-R3*sind(phi)-R4*cosd(o)-R_1y)^2)

(Lmono)^2 == ((((R3*cosd(phi))-(R4*sind(o))-((Ip*cosd(eta o))) (R1x))^2)   (((R3*sind(phi)) (R4*cosd(o))-((Ip*sind(eta o))) (R1y))^2))];

[phi ,o]=vpasolve(eqns,[phi o]);

I wrote this in python :

import math as m

Lb = 145.0
Lm = 320.0
IPS= 16.03;
Pivot= -20.0;
R1x=(-46.72)-((Pivot)*m.sin(IPS))
R1y=((180.26)-((Pivot)*m.cos(IPS)))
R_1x=(-43.52)-((Pivot)*m.sin(IPS))
R_1y=((-151.37)-((Pivot)*m.cos(IPS)))
R4=60.0
R3=203.727
Ip=100.1
eta=36.79

import sympy as sym
from sympy import sin, cos
    
sym.init_printing()
phi,o = sym.symbols('phi,o')
f = sym.Eq(((R3*cos(phi)-R4*sin(o) R_1x)**2) ((-R3*sin(phi)-R4*cos(o)-R_1y)**2),Lb**2)
g = sym.Eq(((((R3*cos(phi))-(R4*sin(o))-((Ip*cos(eta o))) (R1x))**2)   (((R3*sin(phi)) (R4*cos(o))-((Ip*sin(eta o))) (R1y))**2)),Lm**2)
print(sym.nonlinsolve([f,g],(phi,o)))

But when I run the code it loads for about 30 seconds (in matlab it takes 1-2 seconds) and then returns this:

runfile('C:/Users/Administrator/.spyder-py3/temp.py', wdir='C:/Users/Administrator/.spyder-py3')
EmptySet

EmptySet?

can someone help me?

CodePudding user response：

You can use SymPy's nsolve with an initial guess:

In [7]: sym.nsolve([f, g], [phi, o], [0, 0])
Out[7]: 
⎡0.228868116194702⎤
⎢                 ⎥
⎣0.345540167046199⎦

CodePudding user response：

Oscar's answer is "correct", but since you are new to Python there are a few things that you'd want to know.

First, in Matlab you are using sind and cosd which require the angle in degrees. On the other hand, the trigonometric functions exposed by the math, numpy, sympy modules require the angle to be in radians. Hence, we need to convert Lb, Lm, ... to radians.

NOTE: since I don't know the kind of problem you are solving, I have applied the m.radians to every number. This is probably wrong: you have to fix it!

Once we are done with that, we can use nsolve to numerically solve the system of equation, by providing an initial guess.

import math as m
import sympy as sym
from sympy import sin, cos, lambdify, nsolve, Add

sym.init_printing()
phi,o = sym.symbols('phi,o')

Lb = m.radians(145.0)
Lm = m.radians(320.0)
IPS= m.radians(16.03)
Pivot= m.radians(-20.0)
R1x=m.radians(-46.72)-((Pivot)*m.sin(IPS))
R1y=m.radians(180.26)-((Pivot)*m.cos(IPS))
R_1x=m.radians(-43.52)-((Pivot)*m.sin(IPS))
R_1y=m.radians(-151.37)-((Pivot)*m.cos(IPS))
R4=m.radians(60.0)
R3=m.radians(203.727)
Ip=m.radians(100.1)
eta=m.radians(36.79)
f = sym.Eq(((R3*cos(phi)-R4*sin(o) R_1x)**2) ((-R3*sin(phi)-R4*cos(o)-R_1y)**2),Lb**2)
g = sym.Eq(((((R3*cos(phi))-(R4*sin(o))-((Ip*cos(eta o))) (R1x))**2)   (((R3*sin(phi)) (R4*cos(o))-((Ip*sin(eta o))) (R1y))**2)),Lm**2)

print(nsolve([f, g], [phi, o], [0, 0]))
# out: Matrix([[0.560675440923978], [-0.0993239452750302]])

Since you are solving an optimization problem and the equations are non-linear, it is likely that there are more than one solution. We can create a contour plot of the two equations: the intersections between the contours represent the solutions:

# convert symbolic expression to numerical functions for evaluation
fn = lambdify([phi, o], f.rewrite(Add))
gn = lambdify([phi, o], g.rewrite(Add))

# plot the 0-level contour of fn and gn: the intersection between
# the curves are the solutions you are looking for
import numpy as np
import matplotlib.pyplot as plt
pp, oo = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j]
fig, ax = plt.subplots()
ax.contour(pp, oo, fn(pp, oo), levels=[0], cmap="Greens_r")
ax.contour(pp, oo, gn(pp, oo), levels=[0], cmap="winter")
ax.set_xlabel("phi [rad]")
ax.set_ylabel("o [rad]")
plt.show()

From this picture you can see that there are two solutions. We can find the other solution by providing a better initial guess to nsolve:

print(nsolve([f, g], [phi, o], [0.5, 3]))
# out: Matrix([[0.451286281041857], [2.54087384334971]])