Solving system of non-linear trigonometric symbolic equations

2 views (last 30 days)
I have created a matlab function to model the rotation of a particular configuration of free rotating boxes. I end up with a system of symbolic trigonometric equations that I must solve however matlab gives strange answers when I increase the number of boxes to greater than one. After running the mechanics equations the system I get is such where M_Sum is the sum of the moments (equilibrium so they must equal zero):
M_SUM = [(cos(theta1) + sin(theta1)/2)*(cos(theta1)/2 + sin(theta1)) + (cos(theta1) - sin(theta1)/2)*(cos(theta2)/2 - cos(theta1)/2 + sin(theta1) + sin(theta2) - 2/5) - (cos(theta1)/2 + sin(theta1))*(cos(theta1) + cos(theta2) + sin(theta1)/2 - sin(theta2)/2 + 11/10) + (cos(theta1)/2 - sin(theta1))*(cos(theta1) - sin(theta1)/2 + 11/10)
(cos(theta2) + sin(theta2)/2)*(cos(theta2)/2 - cos(theta1)/2 + sin(theta1) + sin(theta2) - 2/5) + (cos(theta2) - sin(theta2)/2)*(cos(theta3)/2 - cos(theta2)/2 + sin(theta2) + sin(theta3) + 1/10) - (cos(theta2)/2 + sin(theta2))*(cos(theta2) + cos(theta3) + sin(theta2)/2 - sin(theta3)/2 + 11/10) + (cos(theta2)/2 - sin(theta2))*(cos(theta1) + cos(theta2) + sin(theta1)/2 - sin(theta2)/2 + 11/10)
(cos(theta3) - sin(theta3)/2)*(cos(theta4)/2 - cos(theta3)/2 + sin(theta3) + sin(theta4) + 2/5) + (cos(theta3) + sin(theta3)/2)*(cos(theta3)/2 - cos(theta2)/2 + sin(theta2) + sin(theta3) + 1/10) - (cos(theta3)/2 + sin(theta3))*(cos(theta3) + cos(theta4) + sin(theta3)/2 - sin(theta4)/2 + 11/10) + (cos(theta3)/2 - sin(theta3))*(cos(theta2) + cos(theta3) + sin(theta2)/2 - sin(theta3)/2 + 11/10)
(cos(theta4) + sin(theta4)/2)*(cos(theta4)/2 - cos(theta3)/2 + sin(theta3) + sin(theta4) + 2/5) - (cos(theta4)/2 + sin(theta4))*(cos(theta4) + sin(theta4)/2 + 41/10) - (cos(theta4) - sin(theta4)/2)*(cos(theta4)/2 - sin(theta4) + 1/10) + (cos(theta4)/2 - sin(theta4))*(cos(theta3) + cos(theta4) + sin(theta3)/2 - sin(theta4)/2 + 11/10)]
This is a system of 4 equations and 4 unknowns. When I use the "solve" command, setting M_Sum==0 (0 vector in this case [0;0;0;0]), it takes a very long time and provides me a struct of answers which aren't values but rather contain "root" and "z" in them. How can I get the theta values theta1->4 by setting M_Sum==0? Eventually I hope to expand this to dozens or more equations and unknowns within the limits of computational speed.
Thank You!

Answers (2)

Stephan
Stephan on 25 Jun 2019
syms theta1 theta2 theta3 theta4
format long
M_SUM = [(cos(theta1) + sin(theta1)/2)*(cos(theta1)/2 + sin(theta1)) + (cos(theta1) - sin(theta1)/2)*(cos(theta2)/2 - cos(theta1)/2 + sin(theta1) + sin(theta2) - 2/5) - (cos(theta1)/2 + sin(theta1))*(cos(theta1) + cos(theta2) + sin(theta1)/2 - sin(theta2)/2 + 11/10) + (cos(theta1)/2 - sin(theta1))*(cos(theta1) - sin(theta1)/2 + 11/10);
(cos(theta2) + sin(theta2)/2)*(cos(theta2)/2 - cos(theta1)/2 + sin(theta1) + sin(theta2) - 2/5) + (cos(theta2) - sin(theta2)/2)*(cos(theta3)/2 - cos(theta2)/2 + sin(theta2) + sin(theta3) + 1/10) - (cos(theta2)/2 + sin(theta2))*(cos(theta2) + cos(theta3) + sin(theta2)/2 - sin(theta3)/2 + 11/10) + (cos(theta2)/2 - sin(theta2))*(cos(theta1) + cos(theta2) + sin(theta1)/2 - sin(theta2)/2 + 11/10);
(cos(theta3) - sin(theta3)/2)*(cos(theta4)/2 - cos(theta3)/2 + sin(theta3) + sin(theta4) + 2/5) + (cos(theta3) + sin(theta3)/2)*(cos(theta3)/2 - cos(theta2)/2 + sin(theta2) + sin(theta3) + 1/10) - (cos(theta3)/2 + sin(theta3))*(cos(theta3) + cos(theta4) + sin(theta3)/2 - sin(theta4)/2 + 11/10) + (cos(theta3)/2 - sin(theta3))*(cos(theta2) + cos(theta3) + sin(theta2)/2 - sin(theta3)/2 + 11/10);
(cos(theta4) + sin(theta4)/2)*(cos(theta4)/2 - cos(theta3)/2 + sin(theta3) + sin(theta4) + 2/5) - (cos(theta4)/2 + sin(theta4))*(cos(theta4) + sin(theta4)/2 + 41/10) - (cos(theta4) - sin(theta4)/2)*(cos(theta4)/2 - sin(theta4) + 1/10) + (cos(theta4)/2 - sin(theta4))*(cos(theta3) + cos(theta4) + sin(theta3)/2 - sin(theta4)/2 + 11/10)];
sol = vpasolve(M_SUM==[0;0;0;0]);
theta = double([sol.theta1;sol.theta2;sol.theta3;sol.theta4])
format short
  6 Comments
John D'Errico
John D'Errico on 25 Jun 2019
Solve IS built to handle nonlinear systems. Except that many nonlinear problems have no analytical solution! Or, if they do, the result is often a very nastily complicated mess. That is why you found that solve took so long. Each of solve, vpasolve, and fsolve are "real" solvers. They just work differently.
Noam Kaplan
Noam Kaplan on 26 Jun 2019
Walter, cell2mat(struct2cell(sol)) isn't working. It throws an error "CELL2MAT does not support cell arrays containing cell arrays or objects." However when I do sol.theta1 it gives me a numeric value so I know that value is there.

Sign in to comment.


Walter Roberson
Walter Roberson on 26 Jun 2019
There is a closed form solution to the equations. It involves arctan of ratios of polynomials of degree 43, involving the roots of a polynomial of degree 44. The polynomial of degree 44 has 20 real roots
That is, there are 44 sets of solutions to your system, at least 20 of which are entirely real-valued. MATLAB is returning to you the sets of solutions. All of them are actual solutions. Because you use solve(), MATLAB returns the closest you can get to closed form solutions -- the exact rational terms involved, combined with a placeholder for where the roots of the degree 44 polynomial go. There are no closed form solutions of the degree 44 polynomial in terms of "algebraic numbers", so it cannot just give you 44 really really long equations.
You can tell MATLAB to find all of the roots of the polynomials in numeric forms, and then substitute those into the formulas returned by solve() in order to find numeric approximations of the roots.
You could reduce the number of solutions by putting constraints on the variables, such as requiring that they all be real-valued or that they all be non-negative. It looks to me as if there might be only one solution that is entirely real-valued and entirely positive,
theta1 = 3.050815176, theta2 = 0.08374495950, theta3 = 0.3152774464, theta4 = 3.010009959
When you use vpasolve() you can pass in ranges of values to search over, such as
vpasolve(M_SUM, [theta1, theta2, theta3, theta4], [0 pi; 0 pi; 0 pi; 0 pi])

Products


Release

R2017a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!