Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.
How to correctly find the solution of a system of equations containing trigonometric functions
Show older comments
theta=70;L=1500;h=200;D=1;
eq=[l1+l2==L;
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/l2;
(cosd(theta)*cosd(x1)+sind(theta)*sind(x1))/x==cosd(x1)/l1;
tand(x1)==D/(2*f);tand(x1)==x/y;cosd(theta)==h/(y+f)];
vars=[x1,l1,l2,f,x,y];
[x1,l1,l2,f,x,y]=solve(eq,vars)
The code is as above，the solution show
l1 =
7590603.4479823424160756070250191
l2 =
-7587603.8886555868349951502850662,I think the solution is wrong,but I don't know what is wrong，Can anyone fix this?
Answers (3)
It would be difficult to plot this since there are so many variables, and it is not obvious (to me, at least) what 1 or 2 variables you would want to plot (as a line plot or a surface).
With trogonometric equations, there can be an infinity of solutions (I am not certain about this system), with those solutions having periodic characteristics. If you know what the approximate solutions should be, perhaps using vpasolve and giving it the approximate solutions as initial estimates would give you the result you want.
syms x1 l1 l2 f x y
theta=70;
L=1500;
h=200;
D=1;
eq=[l1+l2==L;
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/l2;
(cosd(theta)*cosd(x1)+sind(theta)*sind(x1))/x==cosd(x1)/l1;
tand(x1)==D/(2*f);
tand(x1)==x/y;
cosd(theta)==h/(y+f)];
vars=[x1,l1,l2,f,x,y];
[x1,l1,l2,f,x,y]=solve(eq,vars)
.
John D'Errico
on 2 Sep 2025
Edited: John D'Errico
on 2 Sep 2025
It was smart of you to look at the results, and question them. Always think about what comes out from a solver. Does the solution make sense?
I would suggest that, by the way, naming variables l1 and l2 is always a bad idea. Lower case l an easily be mistaken or mistyped, and depending on the font, it often becomes difficult to disentagle from the number 1. Just a suggestion.
Anyway, trigonometric equations essentially always have infinitely many solutions, at least if they have any soutions at all. And since this system will not have any symboic solution (solve gave up when you tried, and it passed the buck on to vpasolve.) But vpasolve will try to find A solution vased on a set of random starting values.
syms x1 l1 l2 f x y
theta=70;L=1500;h=200;D=1;
eq=[ l1+l2==L;
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/l2;
(cosd(theta)*cosd(x1)+sind(theta)*sind(x1))/x==cosd(x1)/l1;
tand(x1)==D/(2*f);
tand(x1)==x/y;
cosd(theta)==h/(y+f)];
vars=[x1,l1,l2,f,x,y];
So it looks like you indeed have 6 equations, and 6 variables. But a system of 6 equations and 6 unknowns becomes quite difficult to visualize what is happening.
[x1,l1,l2,f,x,y]=vpasolve(eq,vars)
So, is this solution correct? Well, maybe not useful or meaningful. It looks like vpasolve is converging to a somewhat degenerate solution. For example, we have
cosd(theta)==h/(y+f)
but if you look at the values for y and f, they are both huge, but essentially almost the same number, except for opposite signs. That means the denominator in that expression will be the result of a massive subtractive cancellation. The denominator in that expression is based on the least significant digits in f and y. And that is a massive red flag to me.
It tends to suggest to me that possibly your system may not be well-posed. Or it may simply be incorrect. Or it may be it is based on some assumptions which are ony approximately correct, so not fully valid. The problem is, we have absolutely no idea where the equations came from, what the variables represent.
If I were you, I would first return to your equations. Are they correct? Did you mistype something? Is there something strange about what you have done? For example, look at your second equation.
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/l2;
On the left hand side, we have an interesting form. Do you recognize the form for the cosine of the sum of two angles?
In there, you will see the identity
cos(alpha +/- beta) = cos(alpha)*cos(bets) -/+ sin(alpha)*sin(beta)
Which suggests you can rewrite equation 2 as:
cosd(theta + x1)/x = cosd(x1)/l2
In my eyes, this is simpler, but I'm not sure it helps.
Anyway, can we fix it? No. That is to a large extent, because we are given no clue as to what the equations mean, what they model, and what the variables represent.
David Goodmanson
on 6 Sep 2025
Edited: David Goodmanson
on 6 Sep 2025
Hello yan,
This is one of those situations where a multivariable solver is best avoided if possible. With enough variable substitution and elimination you can reduce the system to finding the zeros of a function of one variable. I did the algebra by hand which is increasingly falling out of favor but has the advantage of working. I used
z in place of x1
and in your fourth line of equations, absorbed the 2 in D/2 into a new D
(so if old D = 1, then new D = 1/2)
The one-variable function R(z) is defined below, and a plot of R vs z looks like
It repeats every 180 degrees and zooming in showed zeros near 13+ and 151+ degrees which were pinned down using fzero. This produces the the following two solutions:
z12 = 13.679 151.39 % roots in degrees
f = 2.0544 -0.9168
y = 582.71 585.68
x = 141.82 -319.41
i2 = 1251.5 -373.8
i1 = 248.49 1873.8
Details: with trig identities, the second and third equations are
cosd(th+z)/x = cosd(z)/i2
cosd(th-z)/x = cosd(z)/i1
and defining
u(z) = cosd(th+z)/cosd(z)
w(z) = cosd(th-z)/cosd(z)
now the second and third equations are
u(z) = x/i2
w(z) = x/i1
the next two are
tand(z) = D/f = x/y
the next one is
y+f = h/cosd(th) = k defining a new constant k
and we still have
i1+i2 = L
Clearly the z variable is the one not to eliminate, and after a few steps the resulting function of z is
R(z) = (tand(z)*k-D)*(1/u(z) + 1/w(z)) - L
The code is
% plot
z = (0:360)+.2;
R = fun(z);
fig(1)
plot(z,R)
xlabel('z')
ylabel('R')
% find roots
th = 70;
L = 1500;
h = 200;
D = 1/2;
k = h/cosd(th);
z1 = fzero(@fun,[13 14])
z2 = fzero(@fun,[151 152])
format short g
z12 = [z1 z2]
% back substitute
f = D./tand(z12)
y = k-f
x = y*D./f
u = cosd(th+z12)./cosd(z12);
w = cosd(th-z12)./cosd(z12);
i2 = x./u
i1 = x./w
format
% check, all of these are small
cosd(th+z12)./x - cosd(z12)./i2
cosd(th-z12)./x - cosd(z12)./i1
tand(z12) - D./f
tand(z12) - x./y
cosd(th) -h./(y+f)
i1+i2 - L
function R = fun(z)
th = 70;
L = 1500;
h = 200;
D = 1/2;
k = h/cosd(th);
u = cosd(th+z)./cosd(z);
w = cosd(th-z)./cosd(z);
t = tand(z);
R = (t*k-D).*(1./u + 1./w) - L;
6 Comments
Hi David,
Amazing insight.
Inspired by your solution, I decided to try it using the tools at hand to let the computer do the algebra, which might not be as fun, but mitigates possibility of errors in all of the manipulations.
syms x1 I1 I2 f x y real
theta=70;L=1500;h=200;D=1;
The equations are easier to see if use symbolic contants for now
syms theta L h real
%{
eqn=[ I1+I2==L;
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/I2;
(cosd(theta)*cosd(x1)+sind(theta)*sind(x1))/x==cosd(x1)/I1;
tand(x1)==D/(2*f);
tand(x1)==x/y;
cosd(theta)==h/(y+f)];
%}
Define the equations using rad to so we don't have pi/180 clutter, convert to deg at the end
eqn=[ I1+I2==L;
(cos(theta)*cos(x1)-sin(theta)*sin(x1))/x==cos(x1)/I2;
(cos(theta)*cos(x1)+sin(theta)*sin(x1))/x==cos(x1)/I1;
tan(x1)==D/(2*f);
tan(x1)==x/y;
cos(theta)==h/(y+f)];
eqn
Anything in a denominator cannot be zero
assumeAlso([x,y,I2,I1,f+y] ~= 0)
Simplify the equations
eqn = simplify(eqn)
Expressions for I2/x and I1/x
eqn1 = isolate(eqn(2),I2)
eqn1 = lhs(eqn1)/x == rhs(eqn1/x)
eqn2 = isolate(eqn(3),I1);
eqn2 = lhs(eqn2)/x == rhs(eqn2/x)
Combine eqn2 and eqn1
eqn3 = eqn2 + eqn1
eqn3 = simplify(lhs(eqn3)) == simplify(rhs(eqn3))
And sub in eqn(1) to eliminate I1 + I2
eqn3 = subs(eqn3,lhs(eqn(1)),rhs(eqn(1)))
Combine eqn(4) and eqn(5) to get an expression for tan(x1)*(f+y)
eqn4 = eqn(4) - eqn(5)
eqn4 = isolate(eqn4,tan(x1))
eqn4 = isolate(eqn4,x + 1/2)
Use eqn(6) to get an expression for f + y in terms of the knowns
eqn5 = isolate(eqn(6),f+y)
and sub that back into eqn4 to get an expression in terms of x and x1
eqn6 = subs(eqn4,f+y,rhs(eqn5))
Reciprocal of eq3 (shouldn't be necessary, but looked a bit cleaner to me)
eqn7 = 1/eqn3
and sub in the expression for x from eqn6
eqn8 = subs(eqn7,x,rhs(isolate(eqn6,x)))
Now we have an expression for one unknown, x1
Sub in the values for the constants
theta = 70*sym(pi)/180; % rad
L = 1500;
h = 200;
eqn8 = subs(eqn8)
solve actually returns a closed form expression for four solutions, but it's kind of nasty. I didn't try to do anything with it.
x1 = solve(eqn8)
Travel to VPA land
x1 = vpa(x1)
Eliminate the small imaginary part
x1 = real(x1)
Backsolve eqn7 for x based on x1
x = vpa(rhs(subs(isolate(eqn7,x))))
Backsolve eqn1 for I2
I2 = vpa(subs(rhs(isolate(eqn1,I2))))
Backsolve eqn(1) for I1
I1 = subs(rhs(isolate(eqn(1),I1)))
Backsolve eqn(4) and eqn(5) for f and y respectively
f = subs(rhs(isolate(eqn(4),f)))
y = subs(rhs(isolate(eqn(5),y)))
Put all four solutions into an array for easier visualization. We see that solutions 2 and 3, not surprisingly, match your result.
%{
z12 = 13.679 151.39 % roots
f = 2.0544 -0.9168
y = 582.71 585.68
x = 141.82 -319.41
b = 1251.5 -373.8
a = 248.49 1873.8
%}
vpa([x1*180/vpa(pi),f,y,x,I2,I1].',5)
David Goodmanson
on 6 Sep 2025
Edited: David Goodmanson
on 7 Sep 2025
Hi Paul, it's interesting to see how it can be accomplished in symbolics and how useful the isolate function is. Of the solutions, the first is the same as the third and the fourth is the same as the second, with 180 subtracted from the angle in both cases so they are basically duplicate solutions, which is good.
Your observation that there are closed form solutions got me interested to go back for another try, and I realized that the sum-of-angle identities were actually not that helpful. With
c = cos(z) s = sin(z) t = tan(z) % radians
cbar = cos(th) sbar = sin(th) % th = (pi/180)*70
then the expressions I had
u(z) = cos(th+z)/cos(z)
w(z) = cos(th-z)/cos(z)
are actually better as
u(z) = cbar - sbar*t
w(z) = cbar + sbar*t
The function R(z) is (k is as defined in the code I used before)
R = (t*k-D)*(1/u + 1/w) - L
= (t*k-D)*(1/(cbar-sbar*t) + 1/(cbar+sbar*t)) - L;
and from this you can obtain a simple quadratic equation in t:
A = L*sbar^2;
B = 2*cbar*k;
C = -L*cbar^2-2*cbar*D;
t0 = roots([A B C])
angles = (180/pi)*atan(t0)
t0 = -0.5454 0.2434
angles = -28.6068 13.6787
same as before for solution1 ( = solution3 - 180) and solution2.
Paul
on 6 Sep 2025
Kudos to you.
I tried to get a closed-form solution by manipulating eqn8 and almost got it to a quadratic with cos(x1) as the unkown, but I couldn't get rid of a pesky cos(x1)*sin(x1) term. Maybe I'll get back to it.
David Goodmanson
on 7 Sep 2025
The sin*cos term is kind of annoying but I think if you isolate that term on one side of the equation and square the equation you can maybe end up with a quadratic in cos^2 which would be another solution method.
Why didn't I think of that ?! (But there is a complication, see below)
syms x1 I1 I2 f x y real
theta=70;L=1500;h=200;D=1;
syms theta L h real
%{
eqn=[ I1+I2==L;
(cosd(theta)*cosd(x1)-sind(theta)*sind(x1))/x==cosd(x1)/I2;
(cosd(theta)*cosd(x1)+sind(theta)*sind(x1))/x==cosd(x1)/I1;
tand(x1)==D/(2*f);
tand(x1)==x/y;
cosd(theta)==h/(y+f)];
%}
eqn=[ I1+I2==L;
(cos(theta)*cos(x1)-sin(theta)*sin(x1))/x==cos(x1)/I2;
(cos(theta)*cos(x1)+sin(theta)*sin(x1))/x==cos(x1)/I1;
tan(x1)==D/(2*f);
tan(x1)==x/y;
cos(theta)==h/(y+f)];
assumeAlso([x,y,I2,I1,f+y] ~= 0);
eqn = simplify(eqn)
eqn1 = isolate(eqn(2),I2);
eqn1 = lhs(eqn1)/x == rhs(eqn1/x);
eqn2 = isolate(eqn(3),I1);
eqn2 = lhs(eqn2)/x == rhs(eqn2/x);
eqn3 = eqn2 + eqn1;
eqn3 = simplify(lhs(eqn3)) == simplify(rhs(eqn3));
eqn3 = subs(eqn3,lhs(eqn(1)),rhs(eqn(1)));
eqn4 = eqn(4) - eqn(5);
eqn4 = isolate(eqn4,tan(x1));
eqn4 = isolate(eqn4,x + 1/2);
eqn5 = isolate(eqn(6),f+y);
eqn6 = subs(eqn4,f+y,rhs(eqn5));
eqn7 = 1/eqn3;
eqn8 = subs(eqn7,x,rhs(isolate(eqn6,x)));
eqn8 = lhs(eqn8) - rhs(eqn8) == 0;
[num,den] = numden(lhs(eqn8));
Here is the pesky sin(x1)*cos(x1)
eqn9 = subs(num,tan(x1),rewrite(tan(x1),'sincos')) == 0
Following your advice
eqn9 = isolate(eqn9,sin(x1)*cos(x1))
eqn9 = eqn9^2
eqn9 = subs(eqn9,sin(x1)^2,rewrite(sin(x1)^2,'cos'))
Now we get a closed form solution, which isn't too bad
sol = solve(eqn9,x1,'ReturnConditions',true)
sol.x1
Plug in the specific values from the problem statement
theta=70*sym(pi)/180;L=sym(1500);h=sym(200);D=sym(1);
subs(sol.x1)
And make sure the conditions are true for those values
syms k integer
isAlways(subs(sol.conditions))
Let k = 0
x1 = vpa(subs(subs(sol.x1),k,0));
x1*180/vpa(pi)
We see that the third and fourth solutions are the same as two derived previously, but the first and second look like phantom solutions that do not solve the original equations that, I assume, came to be from the squaring operation above.
x = vpa(rhs(subs(isolate(eqn7,x))));
I2 = vpa(subs(rhs(isolate(eqn1,I2))));
I1 = subs(rhs(isolate(eqn(1),I1)));
f = subs(rhs(isolate(eqn(4),f)));
y = subs(rhs(isolate(eqn(5),y)));
for ii = 1:numel(eqn)
check(ii,:) = (vpa(subs(lhs(eqn(ii))-rhs(eqn(ii))))).';
end
vpa(check,5)
David Goodmanson
on 8 Sep 2025
Edited: David Goodmanson
on 8 Sep 2025
Hi Paul, good to see that it works, and there's not much complication really. One nice feature is that the (sigma + k*pi) form shows that the solutions are periodic with period pi. Taking -pi/2 < z <= pi/2 as the fundamental interval, tan has a one-element inverse image z in that interval, while cos^2 has a two-element inverse image, z and pi-z, one of which is bad. It's interesting that if you take one of the spurious z values and do all the back substitutions to come up with f, x, y, i1, i2, the only one of the six equations that doesn't work is i1+i2 = L.
Categories
Find more on Numeric Solvers in Help Center and File Exchange
on 8 Sep 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
