How to correctly find the solution of a system of equations containing trigonometric functions

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?

Answer 1

[Star Strider]

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)

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

x1 = 

517.32017720707302681674046481624

l1 = 

7590603.4479823424160756070250191

l2 = 

f = 

x = 

●

y = 

2.773576752484404417389800180833e+36

.

Answer 2

[John D'Errico]

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)

x1 = 

517.32017720707302681674046481624

l1 = 

7590603.4479823424160756070250191

l2 = 

f = 

x = 

y = 

2.773576752484404417389800180833e+36

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?

https://en.wikipedia.org/wiki/List_of_trigonometric_identities

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.

Answer 3

[David Goodmanson]

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;

[Paul]

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

eqn = 

Anything in a denominator cannot be zero

assumeAlso([x,y,I2,I1,f+y] ~= 0)

Simplify the equations

eqn = simplify(eqn)

eqn = 

Expressions for I2/x and I1/x

eqn1 = isolate(eqn(2),I2)

eqn1 = 

eqn1 = lhs(eqn1)/x == rhs(eqn1/x)

eqn1 = 

eqn2 = isolate(eqn(3),I1);
eqn2 = lhs(eqn2)/x == rhs(eqn2/x)

eqn2 = 

Combine eqn2 and eqn1

eqn3 = eqn2 + eqn1

eqn3 = 

eqn3 = simplify(lhs(eqn3)) == simplify(rhs(eqn3))

eqn3 = 

And sub in eqn(1) to eliminate I1 + I2

eqn3 = subs(eqn3,lhs(eqn(1)),rhs(eqn(1)))

eqn3 = 

Combine eqn(4) and eqn(5) to get an expression for tan(x1)*(f+y)

eqn4 = eqn(4) - eqn(5)

eqn4 = 

eqn4 = isolate(eqn4,tan(x1))

eqn4 = 

eqn4 = isolate(eqn4,x + 1/2)

eqn4 = 

Use eqn(6) to get an expression for f + y in terms of the knowns

eqn5 = isolate(eqn(6),f+y)

eqn5 = 

and sub that back into eqn4 to get an expression in terms of x and x1

eqn6 = subs(eqn4,f+y,rhs(eqn5))

eqn6 = 

Reciprocal of eq3 (shouldn't be necessary, but looked a bit cleaner to me)

eqn7 = 1/eqn3

eqn7 = 

and sub in the expression for x from eqn6

eqn8 = subs(eqn7,x,rhs(isolate(eqn6,x)))

eqn8 = 

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)

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)

x1 = 

Travel to VPA land

x1 = vpa(x1)

x1 = 

Eliminate the small imaginary part

x1 = real(x1)

x1 = 

Backsolve eqn7 for x based on x1

x = vpa(rhs(subs(isolate(eqn7,x))))

x = 

Backsolve eqn1 for I2

I2 = vpa(subs(rhs(isolate(eqn1,I2))))

I2 = 

Backsolve eqn(1) for I1

I1 = subs(rhs(isolate(eqn(1),I1)))

I1 = 

Backsolve eqn(4) and eqn(5) for f and y respectively

f = subs(rhs(isolate(eqn(4),f)))

f = 

y = subs(rhs(isolate(eqn(5),y)))

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)

ans = 

[David Goodmanson]

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]

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]

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.

[Paul]

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)

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

eqn9 = 

Following your advice

eqn9 = isolate(eqn9,sin(x1)*cos(x1))

eqn9 = 

eqn9 = eqn9^2

eqn9 = 

eqn9 = subs(eqn9,sin(x1)^2,rewrite(sin(x1)^2,'cos'))

eqn9 = 

Now we get a closed form solution, which isn't too bad

sol = solve(eqn9,x1,'ReturnConditions',true)

sol = struct with fields:
            x1: [4×1 sym]
    parameters: k
    conditions: [4×1 sym]

sol.x1

ans = 

●

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)

ans = 

And make sure the conditions are true for those values

syms k integer
isAlways(subs(sol.conditions))

ans = 4×1 logical array
   1
   1
   1
   1

Let k = 0

x1 = vpa(subs(subs(sol.x1),k,0));
x1*180/vpa(pi)

ans = 

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)

ans = 

[David Goodmanson]

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.