How to use syms to solve system of equations?

I am trying to use the syms function to solve for the angles in a 4-bar mechanism. There is an input angle that is "known" (theta2), and two resulting output angles (theta3 and theta4). I would like to be able to solve the system of equations for theta3 and theta4 so that I can run a loop to determine what they are based on theta2 from 0 to 2pi radians. It works already if theta2 is assigned a value (ex. 1 rad) and then solving for the numeric values of theta3 and theta4 using syms, however, this means that syms has to run for every value of theta2 I would like to use. This ends up being slow; it would be much faster if I could solve for theta3 and theta4 in terms of themselves and theta2, so that I can just plug a value of theta2 into the solved equations and spit out an answer. I apologize if this is difficult to follow; below is the script that I am working with.

clc
clear all
close all

ao = 3; %Defining given lengths of members
ab = 6;
bc = 8;
co = 10;

syms theta3 theta4 theta2  %Start the equation solve

assume(theta3>0 & theta3<pi); %Defining boundries for theta3 and theta4
assume(theta4>0 & theta4<pi);

ao_ = ao*[cos(theta2) sin(theta2) 0]; %Use the magnitude and direction of each member to define it as a vector
ab_ = ab*[cos(theta3) sin(theta3) 0];
bc_ = bc*[cos(theta4) sin(theta4) 0];
co_ = [co 0 0];

eq1 = ao_(1) + ab_(1) == bc_(1) + co_(1); %Set up two equations to solve for both theta3 and theta4
eq2 = ao_(2) + ab_(2) == bc_(2) + co_(2);

solution = solve(eq1, eq2, theta3, theta4); %Solving the system of equations for theta 3 and theta 4

Warning: Solutions are parameterized by the symbols: [z1, z2], z2. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

solution = [solution.theta3 solution.theta4];

[Walter Roberson]

theta3 = -2*atan(((105*(32*tan(theta2/2) - 5^(1/2)*((11*tan(theta2/2)^2 + 3)*(9*tan(theta2/2)^2 + 49))^(1/2)))/(11*tan(theta2/2)^2 + 35) - 24*tan(theta2/2) + (33*tan(theta2/2)^2*(32*tan(theta2/2) - 5^(1/2)*((11*tan(theta2/2)^2 + 3)*(9*tan(theta2/2)^2 + 49))^(1/2)))/(11*tan(theta2/2)^2 + 35))/(297*tan(theta2/2)^2 + 105))
theta4 = -2*atan((3*(32*tan(theta2/2) - 5^(1/2)*((11*tan(theta2/2)^2 + 3)*(9*tan(theta2/2)^2 + 49))^(1/2)))/(11*tan(theta2/2)^2 + 35))

[Nicholas Palmiero]

Thank you, but I am looking to learn how to go about getting the result, not really the result itself. Can you tell me how you got that output?

[Walter Roberson]

I used solve with returnconditions. The conditions were that a parameter had to be in a particular range and had to match two computations. I used children() to extract the two computations and solve() without conditions, which gave me two sets of formulas. I then substituted a couple of numeric values for theta2 and vpa() and checked which of the pair of formulas satisfied the range conditions.

As I did not check in detail, it might be the case that somewhere along the line the other formula pair might be needed.

Answer 1

[David Goodmanson]

Hi NIcholas,

This can be done without using syms. Syms of course can be slow, and the result in this case is a typically ponderous syms expression that is hard to interpret. The alternate solution below, which is numerical, has a direct interpretation in terms of the geometry. You do have to vary the values of theta2 by means of a for loop, but 100,000 times through the loop takes less than three seconds. With just one value for th2, the code looks like

ao = 3; %Defining given lengths of members
ab = 6;
bc = 8;
co = 10;
th2 = pi/4;                            % for example
ao_ = ao*[cos(th2) sin(th2) 0];        % Use the magnitude and direction of each member ... 
co_ = [co 0 0];

F = (ao_(1)-co_(1));
G = (ao_(2)-co_(2));
A = ab;
B = bc;
p = conv([A, F+i*G],[F-i*G, A]) - [0 B^2 0];
z3 = roots(p);
z3 = z3(imag(z3) >=0);             % same as th3 >= 0
th3 = angle(z3);
z4 = (F+i*G + A*z3)/B; 
th4 = angle(z4);

% check, should be small
ao*cos(th2) + ab*cos(th3) - ( bc*cos(th4) + co_(1) )
ao*sin(th2) + ab*sin(th3) - ( bc*sin(th4) + co_(2) )

As to a derivation, the starting point is

ao*cos(th2) - co_(1) + ab*cos(th3) = bc*cos(th4) 
ao*sin(th2) - co_(2) + ab*sin(th3) = bc*sin(th4)

Where co_ is taken over the the left hand side. You have the th2 dependence in ao_, so define F,G,A,B as in the code and the equations are

F + A*cos(th3) = B*cos(th4)
G + A*sin(th3) = B*sin(th4)

Take the sum of the first equation and i times the second equation to obtain

F+i*G + A*exp(i*th3) = B*exp(i*th4)       (a)

Defining

z3 = exp(i*th3)   z4 = exp(i*th4)

as unit vectors in the complex plane, then

F+i*G + A*z3 = B*z4                       (b)

is a picture in the complex plane of three of the bars in the mechanism. Taking the complex conjugate of (a), then

F-i*G + A*exp(-i*th3) = B*exp(-i*th4) 

multiply this with (a)

(F+i*G + A*exp(i*th3))*(F-i*G + A*exp(-i*th3)) = B^2

This is an equation for th3 only. To solve it use the fact that since z3 is a unit vector,

exp(-i*th3) = 1/z3

plug that in, then multiply the result by z3

(F+i*G + A*z3)*(F-i*G + A/z3) - B^2 = 0 
(F+i*G + A*z3)*((F-i*G)*z3 + A) -B^2*z3 = 0 

This is a quadratic in z3, and using conv for multiplcation of polynomial coefficients,

p = conv([A, F+i*G],[F-i*G, A]) - [0 B^2 0];
z3 = roots(p);

There are two roots, and positive angle for th3 means that z3 has a positive imaginary part. Take that root, use the angle function to get th3 and go back and use (b) to obtain z4. Since the output of the angle function is -pi < theta <pi, the angle function automatically imposes the upper limit on th3.