How to use syms to solve system of equations?

9 views (last 30 days)
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.
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
Walter Roberson on 14 Mar 2022
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.

Sign in to comment.

Answers (1)

David Goodmanson
David Goodmanson on 15 Mar 2022
Edited: David Goodmanson on 17 Mar 2022
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)
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.





Community Treasure Hunt

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

Start Hunting!