'solve' results in incomplete solutions

Hi,

I'm trying to find the equilibrium points of a system described by 4 equations. This is the code:

 clear all;
 clc;
 syms x1 x2 x3 x4 g m1 m2 l1 l2 pi
 eq1 = x2
 eq2 = -((g*(2*m1+m2)*sin(x1)+m2*(g*sin(x1-2*x3)+2*(l2*x4^2+l1*x2^2*cos(x1- x3))*sin(x1-x3)))/(2*l1*(m1+m2-m2*cos(x1-x3)^2)));
 eq3 = x4
 eq4 = (((m1+m2)*(l1*x2^2+g*cos(x1))+l2*m2*x4^2*cos(x1-x3))*sin(x1-x3))/(l2*(m1+m2-m2*cos(x1-x3)^2));
 eq2=subs(eq2, g, 9.81); eq2=subs(eq2, m1, 2); eq2=subs(eq2, m2, 1);   eq2=subs(eq2, l1, 2); eq2=subs(eq2, l2, 1)
 eq4=subs(eq4, g, 9.81); eq4=subs(eq4, m1, 2); eq4=subs(eq4, m2, 1);  eq4=subs(eq4, l1, 2); eq4=subs(eq4, l2, 1)
 res=solve(0==eq1,0==eq2,0==eq3,0==eq4,x1,x2,x3,x4)
 [res.x1 res.x2 res.x3 res.x4]

And the result I get is:

 [  0, 0,  0,  0]
 [  0, 0, pi,  0]
 [ pi, 0,  0,  0]
 [ z2, 0, z3, z4]

I have 2 problems with this result:

1) All the first 3 are valid solutions, however, the 4th row seems gibberish to me. What is the "z" variable?

2) Also, solve skipped one of the obvious solutions (its a double pendulum system) [pi, 0, pi, 0].

In advance, thank you for the replies.

Answer 1

[Star Strider]

The z2, z3, z4 values I believe are because of this (from the solve documentation):

• If the solution of an equation or a system of equations contains parameters, the solver can choose one or more values of the parameters and return the results corresponding to these values. For some equations and systems, the solver returns parameterized solutions without choosing particular values. In this case, the solver also issues a warning indicating the values of parameters in the returned solutions.

Except it didn't issue any warning indicating it was making up parameter symbols when I ran your code. (It has with other code I've run in the past. I have no idea why it didn't with yours.)

Without changing anything else in your code, I added these two lines (after the respective original lines for eq2 and eq4 respectively:

eq2 = simplify(collect(expand(eq2)))

eq4 = simplify(collect(expand(eq4)))

(because the solve function usually does better after such simplifications), and these:

[X1 X2 X3 X4] = solve(0==eq1,0==eq2,0==eq3,0==eq4, x1,x2,x3,x4);
SymMtx = [X1 X2 X3 X4]
VPAMtx = vpa(SymMtx, 6)

and got this solution:

SymMtx =

[     0, 0,                                        0, 0]
[    pi, 0,                                        0, 0]
[     0, 0,                                       pi, 0]
[  pi/2, 0,    2*atan(((2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[  pi/2, 0,  2*atan((- (2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[  pi/2, 0,   -2*atan(((2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[ -pi/2, 0,    2*atan(((2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[ -pi/2, 0,   -2*atan(((2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[  pi/2, 0, -2*atan((- (2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[ -pi/2, 0,  2*atan((- (2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]
[ -pi/2, 0, -2*atan((- (2^(1/2)*2*i)/3 - 1/3)^(1/2)), 0]

that using vpa produced:

VPAMtx =

[       0, 0,                    0, 0]
[ 3.14159, 0,                    0, 0]
[       0, 0,              3.14159, 0]
[  1.5708, 0,   1.5708 + 1.14622*i, 0]
[  1.5708, 0,   1.5708 - 1.14622*i, 0]
[  1.5708, 0, - 1.5708 - 1.14622*i, 0]
[ -1.5708, 0,   1.5708 + 1.14622*i, 0]
[ -1.5708, 0, - 1.5708 - 1.14622*i, 0]
[  1.5708, 0, - 1.5708 + 1.14622*i, 0]
[ -1.5708, 0,   1.5708 - 1.14622*i, 0]
[ -1.5708, 0, - 1.5708 + 1.14622*i, 0]

I'm not certain this was what you were anticipating (or want), but it at least managed to avoid producing a solution with z2, z3, z4, even though it still only produced three real results. If this wasn't what you were expecting, it might be worth the effort to look over your original code and check for errors.

[Daniel]

Thank you for your reply Star Strider.

I appreciate your efforts in getting rid of the [z2,0,z3,z4] answer, I was afraid it might be masking one of the real answers.

My main problem with this is that one of the trivial answers, [pi, 0, pi, 0], is not found as a solution to the system.

It's a double pendulum system, so the trivial equilibriums are [+-pi or 0, 0, +-pi or 0, 0]

You can easily check that [pi, 0, pi, 0] is a solution (zero dynamic in this case) by substituting x1=pi,x2=0,x3=pi,x4=0 in the equations. However, solve fails to find it.

[Star Strider]

I think the solver is doing its best. You may be encountering numeric precision problems, the likely reason for the last eight rows of the matrix. I multiplied row 4 by its complex conjugate, took the square root, and doubled it, getting:

[ 3.1415926535897824578569270670414, 0, 3.8890676724243935777961962291675, 0]

admittedly a bit of a stretch. That's likely as close as the solver is going to get. You might bring this up with TMW Tech Support, since it could be a bug rather than simply an expected problem given the nature of numeric approximation.

ADDITION — When I commented-out the subs statements, I got this symbolic result:

SymMtx =

[     0, 0,                                                        0, 0]
[     0, 0,                                                       pi, 0]
[    pi, 0,                                                        0, 0]
[  pi/2, 0,   2*atan(((m2 - m1 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[  pi/2, 0,  -2*atan(((m2 - m1 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[ -pi/2, 0,   2*atan(((m2 - m1 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[  pi/2, 0,  2*atan((-(m1 - m2 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[ -pi/2, 0,  -2*atan(((m2 - m1 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[ -pi/2, 0,  2*atan((-(m1 - m2 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[  pi/2, 0, -2*atan((-(m1 - m2 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]
[ -pi/2, 0, -2*atan((-(m1 - m2 + 2*(-m1*m2)^(1/2))/(m1 + m2))^(1/2)), 0]

This may give you some clue as to what the problem may be. It seems it's not a numerical approximation problem after all.