vpasolve giving incorrect answers for a system of trigonometric equations. How do I set bounds to the solution?

I'm trying to use vpasolve to solve a system of trigonomic equations. Known variables are Sa1, Sa2, and Sa3. Unknown variables are S, AR, PHI. The code hasn't been giving me the correct results and I can't figure out why. Would assigning boundaries to the solution help? If so, can someone give me an example of how it would work? I assigned fake values for the unknown variables, so I could calculate known variable values to test. %locations of observation points x1 = -9.72; y1 = -5.62; %OBS4 x2 = -4.95; y2 = -5.60; %OBS5 x3 = 6.63; y3 = -3.76; %OBS6 %assign fake values to unknown variables S = 3.4e-5; AR = 10; PHI = 36; %interim variables W1 = ((-x1*sind(PHI) + y1*cosd(PHI)) / (y1*sind(PHI) + x1*cosd(PHI))); W2 = ((-x2*sind(PHI) + y2*cosd(PHI)) / (y2*sind(PHI) + x2*cosd(PHI))); W3 = ((-x3*sind(PHI) + y3*cosd(PHI)) / (y3*sind(PHI) + x3*cosd(PHI))); %calculate known variable values to test Sa1 = S*(sqrt(AR^-1)*(cosd(W1))^2+sqrt(AR)*(sind(W1))^2) Sa2 = S*(sqrt(AR^-1)*(cosd(W2))^2+sqrt(AR)*(sind(W2))^2) Sa3 = S*(sqrt(AR^-1)*(cosd(W3))^2+sqrt(AR)*(sind(W3))^2) The above returns Sa1 = 1.0752e-05, Sa2 = 1.0753e-05, and Sa3 = 1.0894e-05. Having created known variable values with my fake unknown variable values, I'm not able to return the numbers I assigned. The code below is what I'm trying to use in my project. syms S AR PHI %locations of observation points x1 = -9.72; y1 = -5.62; %OBS4 x2 = -4.95; y2 = -5.60; %OBS5 x3 = 6.63; y3 = -3.76; %OBS6 %values for known variables Sa1 = 1.0752e-05; Sa2 = 1.0753e-05; Sa3 = 1.0894e-05; %interim variables W1 = ((-x1*sind(PHI) + y1*cosd(PHI)) / (y1*sind(PHI) + x1*cosd(PHI))); W2 = ((-x2*sind(PHI) + y2*cosd(PHI)) / (y2*sind(PHI) + x2*cosd(PHI))); W3 = ((-x3*sind(PHI) + y3*cosd(PHI)) / (y3*sind(PHI) + x3*cosd(PHI))); %define system of equations and solve EQN = [Sa1 == S*(sqrt(AR^-1)*(cosd(W1))^2+sqrt(AR)*(sind(W1))^2), ... Sa2 == S*(sqrt(AR^-1)*(cosd(W2))^2+sqrt(AR)*(sind(W2))^2), ... Sa3 == S*(sqrt(AR^-1)*(cosd(W3))^2+sqrt(AR)*(sind(W3))^2)]; [S,AR,PHI] = vpasolve(EQN) The above code returns nonsense. Why? I've tried implementing initial parameter guessing to help the solver along, but that doesn't work either, it just returns empty values. This was my attempt: %try alternative definition of equations and set initial parameter guesses EQN1 = Sa1 == S*(sqrt(AR^-1)*(cosd(W1))^2 + sqrt(AR)*(sind(W1))^2); EQN2 = Sa2 == S*(sqrt(AR^-1)*(cosd(W2))^2 + sqrt(AR)*(sind(W2))^2); EQN3 = Sa3 == S*(sqrt(AR^-1)*(cosd(W3))^2 + sqrt(AR)*(sind(W3))^2);

[sol_S, sol_AR, sol_PHI] = vpasolve([EQN1 EQN2 EQN3], [S AR PHI], [1e-5; 10; 50])

Answer 1

[Paul]

Hi Nat,

I didn't try to analyze the equations to be solved, but they appear to be very sensitive to the actual values of Sa1, Sa2, and Sa3. The assignment of values to those variables loses several digits of precision. The code below does not do the assignment, so the LHS of EQN keeps the full precision that can be retained for symbolic constants after converting from double.

%locations of observation points
x1 = -9.72; y1 = -5.62; %OBS4
x2 = -4.95; y2 = -5.60; %OBS5
x3 =  6.63; y3 = -3.76; %OBS6
%assign fake values to unknown variables
S = 3.4e-5; AR = 10; PHI = 36;
%interim variables
W1 = ((-x1*sind(PHI) + y1*cosd(PHI)) / (y1*sind(PHI) + x1*cosd(PHI)));
W2 = ((-x2*sind(PHI) + y2*cosd(PHI)) / (y2*sind(PHI) + x2*cosd(PHI)));
W3 = ((-x3*sind(PHI) + y3*cosd(PHI)) / (y3*sind(PHI) + x3*cosd(PHI)));
%calculate known variable values to test
Sa1 = S*(sqrt(AR^-1)*(cosd(W1))^2+sqrt(AR)*(sind(W1))^2)

Sa1 = 1.0752e-05

Sa2 = S*(sqrt(AR^-1)*(cosd(W2))^2+sqrt(AR)*(sind(W2))^2)

Sa2 = 1.0753e-05

Sa3 = S*(sqrt(AR^-1)*(cosd(W3))^2+sqrt(AR)*(sind(W3))^2)

Sa3 = 1.0894e-05

syms S AR PHI
%locations of observation points
x1 = -9.72; y1 = -5.62; %OBS4
x2 = -4.95; y2 = -5.60; %OBS5
x3 =  6.63; y3 = -3.76; %OBS6
%values for known variables, just keep the full precision values solved above.
%Sa1 = 1.0752e-05;
%Sa2 = 1.0753e-05;
%Sa3 = 1.0894e-05;
%interim variables
W1 = ((-x1*sind(PHI) + y1*cosd(PHI)) / (y1*sind(PHI) + x1*cosd(PHI)));
W2 = ((-x2*sind(PHI) + y2*cosd(PHI)) / (y2*sind(PHI) + x2*cosd(PHI)));
W3 = ((-x3*sind(PHI) + y3*cosd(PHI)) / (y3*sind(PHI) + x3*cosd(PHI)));
%define system of equations and solve
EQN = [Sa1 == S*(sqrt(AR^-1)*(cosd(W1))^2+sqrt(AR)*(sind(W1))^2), ...
       Sa2 == S*(sqrt(AR^-1)*(cosd(W2))^2+sqrt(AR)*(sind(W2))^2), ...
       Sa3 == S*(sqrt(AR^-1)*(cosd(W3))^2+sqrt(AR)*(sind(W3))^2)];

This call to vpasolve solves for the variables in the order define by symvar

symvar(EQN)
sol = vpasolve(EQN)

sol = struct with fields:
     AR: 9.9999999999998449853561092413536
    PHI: -143.99999999999970213585651254036
      S: 0.00003399999999999973577161116267788

What's confusing to me is that specifying the variables in a different order yields a different result. I don't know why that happens.

sol = vpasolve(EQN,[S,AR,PHI])

sol = struct with fields:
      S: 0.0000040816716892196625274775160908436i
     AR: -0.14035454743876028784021737436727
    PHI: -50.73732711215203519786760912156

Specifying the initial search range yields the expected result.

sol = vpasolve(EQN,[S,AR,PHI],[0 1;0 90;0 90])

sol = struct with fields:
      S: 0.00003399999999999973577161116267788
     AR: 9.9999999999998449853561092413536
    PHI: 36.000000000000297864143487459643

[Nat A]

Thank you! I was missing the piece of how to set bounds rather than simply initial guesses. I don't need near as much precision as the results give, but it's of no consequence. This has saved my project from ruin - thanks for the help!

[Paul]

Keep in mind that the bouds specified to vpasolve are bounds on the initial guess for the solution, not necessarily bounds on the solution itself.