Matlab Solve function solutions don't satisfy equation

Question

I have an equation that theoretically has to only have one positive solution and the rest must be negative. So I used matlab solve function but two of the solutions are positive instead of one. Then as expected, by testing the solutions one of them does not satisfy the equation which means it is not the proper results. Here is the code:

n=4
c = 10.3343;
theta = [51.7924 ; 4.3564 ; 0.0103 ; 1.3505];
syms x
eqn = sum(sqrt(x.^2 + theta))  - n.*x == c;
sol = solve(eqn,x);
% find pi
xdouble = double(sol);
xdouble(xdouble<0) = [];

cnew1  = sum(sqrt(xdouble(1).^2 + theta)) - n.*xdouble(1)
cnew2  = sum(sqrt(xdouble(2).^2 + theta)) - n.*xdouble(2)

The results are

cnew1 =

   10.5373


cnew2 =

   10.3343

which clearly second one is the proper answer, but I am confused why matlab generates the first incorrect solution? Any help much appreciated.

Answer 1

It could happens that solve produce spurious solutions (a.k.a incorrect solutions). Check the wikipedia article about spurious solutions.

Why ?

Because solve try to solve your problem analytically. During this process matlab can transform your initial equation to something easier to solve.

For example if we have the following equation:

x+1=0

It is obvious that the solution is x = -1.

But matlab could do something like:

x   = -1
x^2 = (-1)^2 
x^2 = 1

And now we have two solutions: -1 and 1. 1 is a spurious solution.

So you have 2 choices to avoid this problem:

1. Check that the solutions returned by solve are correct.
2. Solve the equation with vpasolve, vpasolve use numerical analysis. So, strictly speaking,vpasolve do not solve your equation but instead determine by iteration the possible solutions.

So:

sol = vpasolve(eqn,x);

will do the trick

Answer 2

1.- MATLAB does not use sum() in this context

In the context of this question MATLAB disregards the sum() operator around sqrt because for sum to be used it must be on a known and finite array of values.

Since an equation is to be solved with solve, and no actual range is defined, because it's command solve the one to find roots, if any root to be found, sum cannot operate without the range defined by a concise array of discrete elements to be summed.

And there's no sampling step either, because for command solve to be accurate solve it has to dive down enough decimals therefore an array step is not specified when defining equations for solve.

So command sum ignored, and in fact one can remove it from the equations to solve.

n=4
c = 10.3343;
th = [51.7924 ; 4.3564 ; 0.0103 ; 1.3505];

figure
ax=gca
hold on

for k=1:1:numel(th)
    f1=@(x)  sum(sqrt(x.^2 + th(k)))  - n.*x - c;
    fplot(ax,f1,[-10 10])
end

plot(ax,[-10 10],[0 0],'Color','k')
grid on

for k=1:1:numel(th)
    f1=@(x)  (sqrt(x.^2 + th(k)))  - n.*x - c;
    fplot(ax,f1,[-10 10])
end

With or without sum the obtained roots are identical :

syms x

xr1=[];

for k=1:1:numel(th)
    
    eqn = sum(sqrt(x.^2 + th(k)))  - n.*x == c;
    
    sol = solve(eqn,x);
    % find pi
    xd = double(sol);
    xr1=[xr1;xd];
    
end

xr1

roots with sum in the equations

xr =
  -0.774025412748041
  -1.881122895419566
  -2.066361839624184
  -2.004354732448066

And without

xr1=[];

for k=1:1:numel(th)
    
    eqn =(sqrt(x.^2 + th(k)))  - n.*x == c;

    sol = solve(eqn,x);
    xd = double(sol);
    xr1=[xr1;xd];

end

xr1

xr =
  -0.774025412748041
  -1.881122895419566
  -2.066361839624184
  -2.004354732448066

Observe that both the shape of the graphs are the same, and the roots are also exactly the same, with or without sum around sqrt

2.- There are 4 equations, not one

Yes one may call an array of equations thrown at solve as just 'an equation' or a 'a vector equation'.

But because the question starts with 'I have an equation' it's worth pointing out that actually there are 4 curves because there are 4 values of theta, I use th, therefore 4 functions, 4 equations, not one, have been defined.

3.- sqrt may have 2 roots

Ok when the parabola only touches the x axis on one point there's only 1 real root, and when the parabola is all above the x axis there are not real roots.

BUT when there are real roots, most likely, they come in pairs.

As poised the equations eqn you are implicitly asking for ANY value satisfying sqrt and that implies that -sqrt should also be considered

for k=1:1:numel(th)
    f1=@(x)  -(sqrt(x.^2 + th(k)))  - n.*x - c;
    fplot(ax,f1,[-10 10])
end

xr2=[];

for k=1:1:numel(th)
    
    eqn = -(sqrt(x.^2 + th(k)))  - n.*x == c;
    
    sol = solve(eqn,x);
    % find pi
    xd = double(sol);
    xr2=[xr2;xd];
    
end

xr2

xr =
  -4.737601253918626
  -3.630503771247101
  -3.445264827042483
  -3.507271934218601

4.- 2nd part of the equations to solve

It seems the cnew expression is also part of the equations, I use cn1 cn2.

Now the resulting cn1 matches the only expected result of 10.3343

cn1=[];
for k=1:1:numel(th)
    cn1  =[cn1 ; (sqrt(xr1(k).^2 + th(k))) - n.*xr1(k)];
end
cn1

cn1 =
  10.334300000000001
  10.334300000000001
  10.334300000000001
  10.334300000000001



cn2=[];
for k=1:1:numel(th)
    cn2  =[cn2 ; (sqrt(xr2(k).^2 + th(k))) - n.*xr2(k)];
end
cn2

cn2 =
  27.566510031349011
  18.709730169976805
  17.227818616339867
  17.723875473748812