strange solution of an equation by Symbolic Toolbox

I come across to an equation I was not able to solve analytically (an advanced math exercise from last year of secondary school of my son). I applied Symbolic Toolbox. It yields two solutions. The 1st one (purely real) fits the equation when put into it. But the 2nd one (purely complex) doesn't ?!? I'm not an expert in Sym. Toolbox, but I follow the help recommendations. Any ideas, recomendations ... ?

the scrip:

%hadanka
clear all
syms x complex    % variable definition
digits(64);     % just for security ...
a=3*sqrt(x^2-9)+4*sqrt(x^2-16)+5*sqrt(x^2-25)==120/x;    % the equation
disp(a);

R=solve(a,x)     % equation solving ...

R = 

x0=vpa(R);     % evaluation of results (following  Help for 'solve' fcn) 
disp(x0);

% check 1st solution
x1=double(x0(1));
L1=3*sqrt(x1^2-9)+4*sqrt(x1^2-16)+5*sqrt(x1^2-25);
R1=120/x1;
disp(num2str([L1 R1 L1-R1]));

24  24   0

% check 2nd solution
x1=double(x0(2));
L2=3*sqrt(x1^2-9)+4*sqrt(x1^2-16)+5*sqrt(x1^2-25);
R2=120/x1;
disp(num2str([L2 R2 L2-R2]));    % PROBLEM !!!

0+74.8552i                 0+26.0707i                 0+48.7845i

% end

Answer 1

[Paul]

This is interesting. Don't understand why solve appears to be using the 12th root in its solution, looks like it should be the 8th

clear all
syms x complex    % variable definition
digits(64);     % just for security ...
a=3*sqrt(x^2-9)+4*sqrt(x^2-16)+5*sqrt(x^2-25)==120/x;    % the equation
disp(a);

R=solve(a,x);     % equation solving ...
x0=vpa(R);     % evaluation of results (following  Help for 'solve' fcn) 
disp(x0);

% check 1st solution
x1=double(x0(1));
L1=3*sqrt(x1^2-9)+4*sqrt(x1^2-16)+5*sqrt(x1^2-25);
R1=120/x1;
disp(num2str([L1 R1 L1-R1]));

24  24   0

% check 2nd solution
format long e
x1=double(x0(2))

x1 = 
      0.000000000000000e+00 - 4.602871649501934e+00i

L2=3*sqrt(x1^2-9)+4*sqrt(x1^2-16)+5*sqrt(x1^2-25);
R2=120/x1;
disp(num2str([L2 R2 L2-R2]));    % PROBLEM !!!

0+74.8552i                 0+26.0707i                 0+48.7845i

% end

Take a look at the second solution.

R(2)

ans = 

Use the 8th root of the polynomial in sigma1 instead of the 12th (I don't know how root orders the roots of the polynomial)

syms z
R2solnew = subs(R(2),root(z^14 + 150*z^12 + 4375*z^10 - 158700*z^8 - 5920625*z^6 + 76593750*z^4 + 1954625625*z^2 - sym("22500000000"), z, 12),root(z^14 + 150*z^12 + 4375*z^10 - 158700*z^8 - 5920625*z^6 + 76593750*z^4 + 1954625625*z^2 - sym("22500000000"), z, 8));

Sub R2solnew into a and check the vpa answer, it seems to check out.

v = vpa(subs(a,x,R2solnew))

v = 

lhs(v) - rhs(v)

ans = 

Check the answer symbolically, checks out

simplify(subs(a,x,R2solnew))

ans = 

symtrue

Check the answer in double, checks out.

format long e
x1=double(R2solnew)

x1 = 
      0.000000000000000e+00 - 2.129977427463911e+00i

L2=3*sqrt(x1^2-9)+4*sqrt(x1^2-16)+5*sqrt(x1^2-25);
R2=120/x1;
([L2 R2 L2-R2]).'    

ans = 
      0.000000000000000e+00 + 5.633862521392059e+01i
      0.000000000000000e+00 + 5.633862521392059e+01i
      0.000000000000000e+00 + 0.000000000000000e+00i

[Petr Kolar]

nice - thank you, but still not very intuitive way I'd say ...

BTW it seems the second solution can be reached by 'fminsearch' too, if starting value x0=2.0i or so ...

[Paul]

No, not intuitive at all. I just hypothesized that perhaps one of the other roots of sigma1 was the correct answer and checked them all. Shouldn't have to do anything like that. You may want to consider filing a bug report with TMW on this problem and see what they say.

[Petr Kolar]

I reported a bug - I'll let you know if any new.

[Petr Kolar]

My collauge proposed a gfaphical solution of the problem: according his approach there are 3 soutions, two of them complex and conjugated.

See graphical plot of error function of the equation.

x = -1:0.05:6;

y = -4:0.05:4;

[xx, yy] = meshgrid(x, y);

zz = xx + i*yy;

a = 3*zz.*sqrt(zz.^2-9)+4*zz.*sqrt(zz.^2-16)+5*zz.*sqrt(zz.^2-25) - 120;

clf();

surf(x, y, abs(a));

[Torsten]

You can square your equation several times.

Finally you will arrive at a polynomial of degree 16 for which you can calculate its roots and check which of these roots satisfies your original equation.

syms x
a = 3*sqrt(x^2-9)+4*sqrt(x^2-16)+5*sqrt(x^2-25)==120/x

a = 

am = expand(a*x)

am = 

am = am - 5*x*sqrt(x^2-25)

am = 

am = simplify(expand(am^2))

am = 

am = am - 12*x^2

am = 

am = simplify(expand(am^2))

am = 

am = am - 2500*x^4 - x^8

am = 

am = expand(am^2)

am = 

am = rhs(am)-lhs(am)

am = 

r = vpasolve(am==0)       % calculate the roots of the deduced polynomial

r = 

a = rhs(a)-lhs(a);
sol = r(abs(double(subs(a,x,r)))<1e-6)   % check which of the roots satisfies your original equation

sol = 

Answer 2

[Torsten]

A MATLAB error.

A second "solution" is computed that does not solve the equation.

MATLAB needs to square your equation several times in order to arrive at the polynomial of degree 14 of which MATLAB computes a root to get a solution of your equation. During this squaring, solutions for the squared systems arise that often don't satisfy the original equation.

But why MATLAB internally does not check the final output - I don't know.

Answer 3

[Petr Kolar]

I can declare with pleasure, that the problem has been solved under R2023a - two solutions are provided: 5.0 and 0-2.19...i (one purely real, other purely imaginary).

Thanks.