The function ''vpasolve'' cannot find a solution with specified range even if its option 'random' was set to true, however the specified range does include 1 solution.

I want to solve nine variables with nine polynomial equations. Firstly, I tried to use the ''solve'' function to obtain all the symbolic solutions. It took a long time to return only one numeric solution with the warning '' Unable to solve symbolically. Returning a numeric solution using vpasolve''. Then, I used the ''vpasolve'' function to find a solution with a specified range but no solution returned. Finally, I used the ''fsolve'' to verify that the specified range includes at least one solution with a appropriate initial value. It is confusing that I cannot find a solution with the specified range by using the ''vpasolve'' function even if the option ''random'' was set to true.

This is the used code:

% Calculation of known coefficients in polynomial equations.
d1=-8;d2=-1;d3=6;
t1=2.2187;t2=1.2392;t3=1.4223;
v=6.038;a1=v*t1/2;a2=v*t2/2;a3=v*t3/2;
c1=d1/2;c2=d2/2;c3=d3/2;
b1=sqrt(a1^2-c1^2);b2=sqrt(a2^2-c2^2);b3=sqrt(a3^2-c3^2);
% nine polynomial equations with nine variables.
syms x1 x2 x3 y1 y2 y3 xr yr r
eq1=(x1-c1)^2/a1^2+y1^2/b1^2-1;
eq2=(x2-c2)^2/a2^2+y2^2/b2^2-1;
eq3=(x3-c3)^2/a3^2+y3^2/b3^2-1;
eq4=(x1-xr)^2+(y1-yr)^2-r^2;
eq5=(x2-xr)^2+(y2-yr)^2-r^2;
eq6=(x3-xr)^2+(y3-yr)^2-r^2;
eq7=c1^2*x1*y1-a1^2*xr*y1+b1^2*x1*yr+b1^2*c1*y1-b1^2*c1*yr;
eq8=c2^2*x2*y2-a2^2*xr*y2+b2^2*x2*yr+b2^2*c2*y2-b2^2*c2*yr;
eq9=c3^2*x3*y3-a3^2*xr*y3+b3^2*x3*yr+b3^2*c3*y3-b3^2*c3*yr;
% equation solving
S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,r,[0 2],'Random',true);
% check if the specified range includes at least one solution
Solution=[1.282855,-3.30307539,1.4754,-3.1486,2.144025872,-3.010426,2,-4,1]
Error=[];
Error(:,1)=subs(eq1,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,2)=subs(eq2,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,3)=subs(eq3,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,4)=subs(eq4,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,5)=subs(eq5,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,6)=subs(eq6,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,7)=subs(eq7,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,8)=subs(eq8,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);
Error(:,9)=subs(eq9,{x1,y1,x2,y2,x3,y3,xr,yr,r},solution);

Answer 1

[John D'Errico]

d1=-8;d2=-1;d3=6;
t1=2.2187;t2=1.2392;t3=1.4223;
v=6.038;a1=v*t1/2;a2=v*t2/2;a3=v*t3/2;
c1=d1/2;c2=d2/2;c3=d3/2;
b1=sqrt(a1^2-c1^2);b2=sqrt(a2^2-c2^2);b3=sqrt(a3^2-c3^2);
% nine polynomial equations with nine variables.
syms x1 x2 x3 y1 y2 y3 xr yr r
eq1=(x1-c1)^2/a1^2+y1^2/b1^2-1;
eq2=(x2-c2)^2/a2^2+y2^2/b2^2-1;
eq3=(x3-c3)^2/a3^2+y3^2/b3^2-1;
eq4=(x1-xr)^2+(y1-yr)^2-r^2;
eq5=(x2-xr)^2+(y2-yr)^2-r^2;
eq6=(x3-xr)^2+(y3-yr)^2-r^2;
eq7=c1^2*x1*y1-a1^2*xr*y1+b1^2*x1*yr+b1^2*c1*y1-b1^2*c1*yr;
eq8=c2^2*x2*y2-a2^2*xr*y2+b2^2*x2*yr+b2^2*c2*y2-b2^2*c2*yr;
eq9=c3^2*x3*y3-a3^2*xr*y3+b3^2*x3*yr+b3^2*c3*y3-b3^2*c3*yr;

Solution=[1.282855,-3.30307539,1.4754,-3.1486,2.144025872,-3.010426,2,-4,1]
Error=[];
Error(:,1)=subs(eq1,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,2)=subs(eq2,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,3)=subs(eq3,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,4)=subs(eq4,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,5)=subs(eq5,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,6)=subs(eq6,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,7)=subs(eq7,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,8)=subs(eq8,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);
Error(:,9)=subs(eq9,{x1,y1,x2,y2,x3,y3,xr,yr,r},Solution);

Error

So, is it a solution? NO. It may be close.

Did you try using vpasolve with that starting value? Did you even call vpasolve properly? No. This is what you wrote:

S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,r,[0 2],'Random',true);

So you listed 9 equations, then the variable r. Why? You have 9 unknowns. Do you somehow expect vpasolve to know what you want to do there, that maybe r is supposed to be between 0 and 2, but that all of the other unknowns are unbounded?

Instead, when I call it properly, it gives me this:

S = 
  struct with fields:

    x1: 1.60884416055053498270434418479
    y1: 2.93704136297912057296719259309
    x2: 1.68004206260311100880431811852
    y2: 3.01305025481464511049150229986
    x3: 1.66160901002679744430328532934
    y3: 2.91905286147157807275399204926
    xr: 1.64960189495632328287532085846
    yr: 2.97021355118473376313175064265
     r: -0.0525508038401069770271615103627

Starting values were not even needed. How well did it do?

     
subs(eq1,S)
ans =
0.0

subs(eq2,S)
ans =
0.0

subs(eq3,S)
ans =
0.0

subs(eq4,S)
ans =
-1.2397791981328813560607768166e-39

And you can see now subs thinks the result is truly zero. Yes, for that last one, 10^-39 is still zero, as far as vpasolve would care.

Is r between 0 and 2? No. That may mean that no solution exists for r in that interval. We don't know, yet. But is your original solution truly a solution? No.

Looking more carefully at your equations. Is it even relevant that r be positive? OF COURSE NOT. r appears in only 3 of those equations.

vpa(eq4,4)
ans =
(x1 - 1.0*xr)^2 + (y1 - 1.0*yr)^2 - 1.0*r^2

vpa(eq5,4)
ans =
(x2 - 1.0*xr)^2 + (y2 - 1.0*yr)^2 - 1.0*r^2

vpa(eq6,4)
ans =
(x3 - 1.0*xr)^2 + (y3 - 1.0*yr)^2 - 1.0*r^2

And in all 3 equations, you square r.

So -0.05 for r is as good as +0.05.

[润晨 李]

Thanks for your careful and helpful answer. As you answered, a reasonable solution can be obtained when starting value are not given and all the variables are unbounded.

S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9);

warning: Solution may be lost.

S =

struct with fields:

x1: 1.6088441605505349827043441847927

y1: 2.9370413629791205729671925930877

x2: 1.6800420626031110088043181185245

y2: 3.0130502548146451104915022998614

x3: 1.6616090100267974443032853293411

y3: 2.9190528614715780727539920492594

xr: 1.6496018949563232828753208584557

yr: 2.9702135511847337631317506426532

r: -0.052550803840106977027161510362718

Be grateful for reminding me to notice that r can be positive or negative. I do admit that the solution meets my stated needs but it is still not what I ultimately want. So I set the option ''Random'' to true in order to obtain more solutions where r is between 0 and 2.

S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,r,[0 2],'Random',true);

It is very disappointing that no solution returned, even for r=0.052550803840106977027161510362718

S =

Empty sym: 0-by-1

Why do I specify the range only for r? In my opinion, it is not necessary to define the ranges for more unknowns when there is no solution only for the range of one unknow. Does the function ''vpasolve'' work better by defining the solution ranges for all unknows?

I tried to solve with the different starting values after hearing form your answer.

x0=[0;0;0;0;0;0;0;0;0]

S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,[x1 y1 x2 y2 x3 y3 xr yr r],x0);

or

[1.282855;-3.30307539;1.4754;-3.1486;2.144025872;-3.010426;2;-4;1]

S=vpasolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,[x1 y1 x2 y2 x3 y3 xr yr r],x0);

No matter how close the starting value was close to the expected solution, the unwanted solution always seemed to return.

S =

struct with fields:

x1: 1.6088441605505349827043441847927

y1: 2.9370413629791205729671925930877

x2: 1.6800420626031110088043181185245

y2: 3.0130502548146451104915022998614

x3: 1.6616090100267974443032853293411

y3: 2.9190528614715780727539920492594

xr: 1.6496018949563232828753208584557

yr: 2.9702135511847337631317506426532

r: -0.052550803840106977027161510362718

The function ''fsolve'' was used to verify the expected solution.

function F=myfun(X)

d1=-8;

d2=-1;

d3=6;

t1=2.2187;

t2=1.2392;

t3=1.4223;

v=6.038;

a1=v*t1/2;

a2=v*t2/2;

a3=v*t3/2;

c1=d1/2;

c2=d2/2;

c3=d3/2;

b1=sqrt(a1^2-c1^2);

b2=sqrt(a2^2-c2^2);

b3=sqrt(a3^2-c3^2);

x1=X(1);

y1=X(2);

x2=X(3);

y2=X(4);

x3=X(5);

y3=X(6);

xr=X(7);

yr=X(8);

r=X(9);

F1=(x1-c1)^2/a1^2+y1^2/b1^2-1;

F2=(x2-c2)^2/a2^2+y2^2/b2^2-1;

F3=(x3-c3)^2/a3^2+y3^2/b3^2-1;

F4=(x1-xr)^2+(y1-yr)^2-r^2;

F5=(x2-xr)^2+(y2-yr)^2-r^2;

F6=(x3-xr)^2+(y3-yr)^2-r^2;

F7=c1^2*x1*y1-a1^2*xr*y1+b1^2*x1*yr+b1^2*c1*y1-b1^2*c1*yr;

F8=c2^2*x2*y2-a2^2*xr*y2+b2^2*x2*yr+b2^2*c2*y2-b2^2*c2*yr;

F9=c3^2*x3*y3-a3^2*xr*y3+b3^2*x3*yr+b3^2*c3*y3-b3^2*c3*yr;

F=[F1;F2;F3;F4;F5;F6;F7;F8;F9];

end

options=optimset('Display','iter','TolFun',1e-16);

x0=[1.282855,-3.30307539,1.4754,-3.1486,2.144025872,-3.010426,2,-4,1];

[x,fval,exitflag]=fsolve('myfun',x0,options);

x =

1.2831 -3.3028 1.4755 -3.1485 2.1437 -3.0104 1.9997 -3.9992 0.9992

The solving error is around 1e-13.

Error=myfun(x);

Error =

1.0e-13 * 0 0 0 -0.0022 0.0011 -0.0011 -0.5684 0.1066 -0.1421

I don't know why the function ''vpasolve'' returns the completely different solution with the same starting values compared to ''fsolve''. I prefer to get as many solutions as possible by specifying the ranges of one or more unknowns rather than using different starting values. It is not easy to select the proper values. Could you give me more good advice?

[John D'Errico]

You say you don't know why fsolve returns a different solution than vpasolve. You need to understand several things.

First, nonlinear problems OFTEN have multiple solutions, each equally as valid as the others. That is certainly true for root finding, which this is.

Do you understand that? The zero a solver finds is not unique, since others will exist. And this is certainly a nonlinear problem. You can talk about uniqueness only for SOME linear problems, and not even all of them.

Next, different algorithms may easily converge to different solutions, when different solutions exist. fsolve and vpasolve are NOT the same algorithm. Worse, a 9 dimensional search space is fairly huge. You can easily hide just about anything in 9 dimensions.

You say that you prefer to get as many solutions as possible. But you don't understand numerical rootfinders. They can NEVER return more than one solution per call. If you start them at a significantly different point, they may find another. Or they may just get stuck, with no root being found.

You have some other questions. But much of what you said seems to worry about r being negative. You do seem to understand that r appears only as the square of r. And that means the sign of r is COMPLETLY, TOTALLY, ABSOLUTELY IRRELEVANT. If you don't like a negative sign there, then take the absolute value of whatever comes out. It will change nothing.

As far as being able to specify a range for only one variable, I saw no explicit comment in the docs for VPASOLVE saying you CAN do that. So I must presume you cannot do that.