How to solve 4 transcendental equations of 4 unknown values?

Hi. I have a problem in solving the equations. I have four equations and four unknowns and I have to find those 4 unknown variables. The equations are as follows:

syms Y2 Y3 l2 l3;
2*Y1*tan(t1_1)+Y2*tan(t2_1)+Y3*tan(t3_1)==0;
2*Y1*tan(t1_3)+Y2*tan(t2_3)+Y3*tan(t3_3)==0;
b1== (t1_1*Y1)+(t2_1*(Y2/2)*((sec(t2_1)^2)/(sec(t1_1)^2)))+(t3_1*(Y3/2)*((sec(t3_1)^2)/(sec(t1_1)^2)));
b3== (t1_3*Y1)+(t2_3*(Y2/2)*((sec(t2_3)^2)/(sec(t1_3)^2)))+(t3_3*(Y3/2)*((sec(t3_3)^2)/(sec(t1_3)^2)));

Here Y1, b1, b3, l1, lambda are known values and the unknown variables are Y2, Y3, l2, l3. The above t1 ,t2 ,t3 are in terms of l2 and l3.

Y1=0.02;
l1 = 0.0105;
lambda =[0.0426  0.0400  0.0405  0.0420];
b1=0.1079;
b3=0.03189;
t2_1=(2*pi/lambda(1))*l2; 
t3_1=(2*pi/lambda(1))*l3; 
t2_3=(2*pi/lambda(3))*l2; 
t3_3=(2*pi/lambda(3))*l3; 
t1_1=(2*pi/lambda(1))*l1;
t1_3=(2*pi/lambda(3))*l1;

I tried using vpasolve and also defining them as functions for plotting contour graphs. But I'm not able to find a solution of all positive values. Please help me with this. Thank you.

Answer 1

[darova]

Use syms and solve

syms Y2 Y3 l2 l3
eq1 = 2*Y1*tan(t1_1)+Y2*tan(t2_1)+Y3*tan(t3_1);
eq2 = 2*Y1*tan(t1_3)+Y2*tan(t2_3)+Y3*tan(t3_3);
eq3 = b1 - (t1_1*Y1)+(t2_1*(Y2/2)*((sec(t2_1)^2)/(sec(t1_1)^2)))+(t3_1*(Y3/2)*((sec(t3_1)^2)/(sec(t1_1)^2)));
eq4 = b3 - (t1_3*Y1)+(t2_3*(Y2/2)*((sec(t2_3)^2)/(sec(t1_3)^2)))+(t3_3*(Y3/2)*((sec(t3_3)^2)/(sec(t1_3)^2)));
solve([eq1 eq2 eq3 eq4])

[LAKKIMSETTI SUNANDA]

Tried this but it's not giving the correct output. I used solve and vpasolve too. But I don't understand why it's unable to show the correct values. Is there any other way to solve these equations?

[Walter Roberson]

syms Y1 l1 b1 b3 
lambda = sym('lambda', [1 3])

lambda = 

syms Y2 Y3 l2 l3
t2_1=(2*pi/lambda(1))*l2; 
t3_1=(2*pi/lambda(1))*l3; 
t2_3=(2*pi/lambda(3))*l2; 
t3_3=(2*pi/lambda(3))*l3; 
t1_1=(2*pi/lambda(1))*l1;
t1_3=(2*pi/lambda(3))*l1;
eq1 = 2*Y1*tan(t1_1)+Y2*tan(t2_1)+Y3*tan(t3_1);
eq2 = 2*Y1*tan(t1_3)+Y2*tan(t2_3)+Y3*tan(t3_3);
eq3 = b1 - (t1_1*Y1)+(t2_1*(Y2/2)*((sec(t2_1)^2)/(sec(t1_1)^2)))+(t3_1*(Y3/2)*((sec(t3_1)^2)/(sec(t1_1)^2)));
eq4 = b3 - (t1_3*Y1)+(t2_3*(Y2/2)*((sec(t2_3)^2)/(sec(t1_3)^2)))+(t3_3*(Y3/2)*((sec(t3_3)^2)/(sec(t1_3)^2)));
eqn = [eq1; eq2; eq3; eq4];
char(eqn)

ans = '[2*Y1*tan((2*l1*pi)/lambda1) + Y2*tan((2*l2*pi)/lambda1) + Y3*tan((2*l3*pi)/lambda1); 2*Y1*tan((2*l1*pi)/lambda3) + Y2*tan((2*l2*pi)/lambda3) + Y3*tan((2*l3*pi)/lambda3); b1 - (2*Y1*l1*pi)/lambda1 + (Y2*l2*pi*cos((2*l1*pi)/lambda1)^2)/(lambda1*cos((2*l2*pi)/lambda1)^2) + (Y3*l3*pi*cos((2*l1*pi)/lambda1)^2)/(lambda1*cos((2*l3*pi)/lambda1)^2); b3 - (2*Y1*l1*pi)/lambda3 + (Y2*l2*pi*cos((2*l1*pi)/lambda3)^2)/(lambda3*cos((2*l2*pi)/lambda3)^2) + (Y3*l3*pi*cos((2*l1*pi)/lambda3)^2)/(lambda3*cos((2*l3*pi)/lambda3)^2)]'

sol1 = solve(eqn(1:2), [Y2,Y3]);
sol1.Y2

ans = 

sol1.Y3

ans = 

eqn2 = subs(eqn(3:4), {Y2,Y3}, {sol1.Y2, sol1.Y3})

eqn2 = 

symvar(eqn2)

ans = 

We have reduced the system to two equations in l2 and l3, given Y1, b1, b3, l1, lambda1, lambda3 . With numeric values of those, there would be the potential for vpasolve() to find a solution. However, you will not be able to find a symbolic solution.

[LAKKIMSETTI SUNANDA]

Thank you but this one is also not giving me the correct values. Can we find the solutions graphically?

[LAKKIMSETTI SUNANDA]

I need a numeric solution. I need to get the numeric values of Y2, Y3, l2 and l3.

[Walter Roberson]

You can fsurf(eqn2) and look for the points that are zero on both surfaces. It is usually a real challenge to interpret.

often it is easier to find the right general area by using

F=sum(eqn2.^2) 
fsurf(F)

This will be 0 anywhere that both equations are 0, and conversely anywhere it is not zero, at least one or both of the two equations are not zero. However it never crosses zero so you can end up spending a long time figuring out whether it touches zero or not.

[Walter Roberson]

You have not given us constants to experiment with, and have not given us the answers you are expecting, so it is difficult for us to test.

[LAKKIMSETTI SUNANDA]

Oh! Sorry. I updated the values. Thank you for your reply.

[darova]

@Walter Roberson he has 4 uknowns

@LAKKIMSETTI SUNANDA I suggest you to express Y1 and use isosurface

[Walter Roberson]

I showed how to reduce the four equations into two equations. Once you have solutions to the two equations you can back-substitute to find the other two values.

Answer 2

[Walter Roberson]

format long g
Q = @(v) sym(v);
 
syms Y2 Y3 l2 l3

Y1 = Q(0.02);
l1 = Q(0.0250);
lambda =[Q(0.1035), Q(0.0999), Q(0.0966)];
b1 = Q('11698461397958468980052614032629')/Q(10)^32;
b3 = Q('12515236662420125800636060478832')/Q(10)^32;

t2_1=(2*pi/lambda(1))*l2; 
t3_1=(2*pi/lambda(1))*l3; 
t2_3=(2*pi/lambda(3))*l2; 
t3_3=(2*pi/lambda(3))*l3; 
t1_1=(2*pi/lambda(1))*l1;
t1_3=(2*pi/lambda(3))*l1;
eq1 = 2*Y1*tan(t1_1)+Y2*tan(t2_1)+Y3*tan(t3_1);
eq2 = 2*Y1*tan(t1_3)+Y2*tan(t2_3)+Y3*tan(t3_3);
eq3 = b1 - (t1_1*Y1)+(t2_1*(Y2/2)*((sec(t2_1)^2)/(sec(t1_1)^2)))+(t3_1*(Y3/2)*((sec(t3_1)^2)/(sec(t1_1)^2)));
eq4 = b3 - (t1_3*Y1)+(t2_3*(Y2/2)*((sec(t2_3)^2)/(sec(t1_3)^2)))+(t3_3*(Y3/2)*((sec(t3_3)^2)/(sec(t1_3)^2)));
eqn = [eq1; eq2; eq3; eq4];
char(eqn)

ans = '[tan((100*pi)/207)/25 + Y2*tan((4000*l2*pi)/207) + Y3*tan((4000*l3*pi)/207); Y2*tan((10000*l2*pi)/483) - tan((233*pi)/483)/25 + Y3*tan((10000*l3*pi)/483); (2000*Y2*l2*pi*cos((100*pi)/207)^2)/(207*cos((4000*l2*pi)/207)^2) - (2*pi)/207 + (2000*Y3*l3*pi*cos((100*pi)/207)^2)/(207*cos((4000*l3*pi)/207)^2) + 11698461397958468980052614032629/100000000000000000000000000000000; (5000*Y2*l2*pi*cos((233*pi)/483)^2)/(483*cos((10000*l2*pi)/483)^2) - (5*pi)/483 + (5000*Y3*l3*pi*cos((233*pi)/483)^2)/(483*cos((10000*l3*pi)/483)^2) + 782202291401257862539753779927/6250000000000000000000000000000]'

sol1 = solve(eqn(1:2), [Y2,Y3]);
sol1.Y2

ans = 

sol1.Y3

ans = 

eqn2 = subs(eqn(3:4), {Y2,Y3}, {sol1.Y2, sol1.Y3})

eqn2 = 

N = 10;
sols = zeros(N,4);
for K = 1 : N
   initial = randn(2,1);  %COLUMN
    sol2 = vpasolve(eqn2, [l2, l3], initial);
    backY2 = subs(sol1.Y2, sol2);
    backY3 = subs(sol1.Y3, sol2);
    sols(K,:) = vpa([backY2, backY3, sol2.l2, sol2.l3]);
end
sols

sols = 10×4

0.126485933403849 -0.167119007948023 -0.220120421257385 0.383807863154484 0.432341841059124 -0.233227313778012 1.13396636173026 0.796286515669898 -2.21198116937819 2.13838491109613 1.69346671322037 1.0702137362349 0.626751243150408 0.493888638447776 -1.66870698078474 1.9101114699078 -0.676916970037412 -0.429495987828577 0.773842328718426 0.22518335023655 0.17078702060627 -0.898235237721649 0.499807867212554 1.24946279610679 -0.872625762873269 -0.405122824760962 0.37514697209635 0.823405507279456 -0.944817212856229 -0.143244176682331 0.474651563810915 1.05004207412011 0.77749208639984 1.12244937333127 1.64897703075677 -2.02387843393713 0.43109071095526 -0.441555710731303 2.07339632941486 0.949422892891187

[LAKKIMSETTI SUNANDA]

Please check the edited code. I tried to plot contours for eq3 and eq4 by substituting Y2,Y3 as you have reduced four equations to two equations. I'm able to see infinite solutions. But I need only positive values of all the variables. Can you tell me how I can see the roots?

[Walter Roberson]

syms Y2 Y3 l2 l3;


Y1=0.02;
l1 = 0.0105;
lambda =[0.0426  0.0400  0.0405  0.0420];
b1=0.1079;
b3=0.03189;
t2_1=(2*pi/lambda(1))*l2; 
t3_1=(2*pi/lambda(1))*l3; 
t2_3=(2*pi/lambda(3))*l2; 
t3_3=(2*pi/lambda(3))*l3; 
t1_1=(2*pi/lambda(1))*l1;
t1_3=(2*pi/lambda(3))*l1;

eqn = [
2*Y1*tan(t1_1)+Y2*tan(t2_1)+Y3*tan(t3_1)==0;
2*Y1*tan(t1_3)+Y2*tan(t2_3)+Y3*tan(t3_3)==0;
b1== (t1_1*Y1)+(t2_1*(Y2/2)*((sec(t2_1)^2)/(sec(t1_1)^2)))+(t3_1*(Y3/2)*((sec(t3_1)^2)/(sec(t1_1)^2)));
b3== (t1_3*Y1)+(t2_3*(Y2/2)*((sec(t2_3)^2)/(sec(t1_3)^2)))+(t3_3*(Y3/2)*((sec(t3_3)^2)/(sec(t1_3)^2)));
];

sol1 = solve(eqn(1:2), [Y2,Y3]);

eqn2 = subs(eqn(3:4), {Y2,Y3}, {sol1.Y2, sol1.Y3});

N = 100;
sols = zeros(N,4);

solcount = 0;
for K = 1 : N
   initial = randn(2,1);  %COLUMN
    sol2 = vpasolve(eqn2, [l2, l3], initial);
    if isempty(sol2) || isempty(sol2.l2); continue; end
    backY2 = subs(sol1.Y2, sol2);
    backY3 = subs(sol1.Y3, sol2);
    thissol = vpa([backY2, backY3, sol2.l2, sol2.l3]);
    if ~all(thissol >= 0); continue; end
    solcount = solcount + 1;    
    sols(solcount,:) = thissol;
    break
end

sols = sols(1:solcount,:);

fprintf('Found %d solutions in %d tries\n', solcount, K);

disp(sols)

All of the candidate solutions I found had at least one negative value; I did not find any solutions that were all positive.

[LAKKIMSETTI SUNANDA]

But it definitely gives positive solutions which are very small like 0.008 to 0.05. I don't understand why its not able to show the solutions of small values. If I try using surf or mesh then I'm able to see those small roots. But its difficult to see each and every solution by clicking at each point. Can you please check if there is any other method to get those roots?

[Alex Sha]

Hi, I don't think there is a solution with all positive values. Be sure the equations you given are all correct.

[Walter Roberson]

I tried over 1300 starting points but none led to all-positive solutions

[LAKKIMSETTI SUNANDA]

Can you try for this input? Can you please check if its giving any positive solution or not?

Y1=0.02;
l1 = 0.0172;
lambda =[0.0712  0.0687  0.0664];
b1=0.1169;
b3=0.1252;
t2_1=(2*pi/lambda(1))*l2; 
t3_1=(2*pi/lambda(1))*l3; 
t2_3=(2*pi/lambda(3))*l2; 
t3_3=(2*pi/lambda(3))*l3; 
t1_1=(2*pi/lambda(1))*l1;
t1_3=(2*pi/lambda(3))*l1;

[Walter Roberson]

>> t4_new
Found 8 solutions in 100 tries
         0.438165660710091         0.480801689969913         0.632205285640874         0.242231767504284
         0.446796293962697         0.087824265522318         0.309580080744266          1.09280300501522
         0.750936013484079          8.92311017431664         0.397348554308343        0.0341199410568209
         0.470778846403452       0.00364711672028139         0.309051631452539          1.90672761491187
          3.23520883575707         0.629358841915416         0.102703543576618         0.327651798541219
         0.509669509078604         0.416105653280171         0.206433666905153         0.703594906946531
         0.499598458592412         0.154187148366108         0.274830616021724          1.20163003959803
        0.0320631935463564         0.331452711890686         0.246724428985903          0.34292876603606

[LAKKIMSETTI SUNANDA]

@Alex Sha Now this is giving us the positive solutions. So I don't think the equations are wrong as I checked those and derived them by hand.

@Walter Roberson We are now getting solutions for a different input. So it should definitely have positive solutions for the above given input too. Those roots will be very close to 0. So if we are able to consider a range of 0 to 0.2, I think we can find them. I'm still working on that and your help will be highly appreciated.

[Walter Roberson]

For the first set of inputs, look at pretty(sol1.Y2)

             / 233 pi \    / 4000 pi l3 \      / 100 pi \    / 10000 pi l3 \
          tan| ------ | tan| ---------- | + tan| ------ | tan| ----------- |
             \   483  /    \     207    /      \   207  /    \     483     /
- ----------------------------------------------------------------------------------
  /    / 4000 pi l2 \    / 10000 pi l3 \      / 4000 pi l3 \    / 10000 pi l2 \ \
  | tan| ---------- | tan| ----------- | - tan| ---------- | tan| ----------- | | 25
  \    \     207    /    \     483     /      \     207    /    \     483     / /

and compare to pretty(sol1.Y3)

           / 233 pi \    / 4000 pi l2 \      / 100 pi \    / 10000 pi l2 \
        tan| ------ | tan| ---------- | + tan| ------ | tan| ----------- |
           \   483  /    \     207    /      \   207  /    \     483     /
----------------------------------------------------------------------------------
/    / 4000 pi l2 \    / 10000 pi l3 \      / 4000 pi l3 \    / 10000 pi l2 \ \
| tan| ---------- | tan| ----------- | - tan| ---------- | tan| ----------- | | 25
\    \     207    /    \     483     /      \     207    /    \     483     / /

The denominators are the same and the numerators have the same form except that one involves l2 and the other involves l3, and there is a negative sign involved.

In order for both of them to be positive, both numerators must have the same sign, and that sign has to be the same as the denominator. Which implies that the sign of the terms with l3 must be negative of the sign of the terms with l2 for the same structure.

This is not impossible, but it is going to take some work to find.