Error using mupadmex Error in MuPAD command: symbolic:s​ym:isAlway​s:LiteralC​ompare|0 < root(z^20 - (555455*z^19)/28224 + (21279533*z^18)/112896 - (524070713*z^17)/451584 + (331

syms r;

eq = (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^2 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^4 * sqrt(((1/2) + 1)^2 - r^2)^2 * r^4 ...
    + (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^8 ...
    - 2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) + 1)^2 - r^2)^2 * sqrt(((1/2) - 1)^2 - r^2)^2 * r^8 ...
    + (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 * r^8;

sol = solve(eq, r);
positive_solutions = sol(sol > 0);

Error using symengine
Unable to prove '0 < root(z^20 - (555455*z^19)/28224 + (21279533*z^18)/112896 - (524070713*z^17)/451584 + (331982393245*z^16)/65028096 - (1096679671049*z^15)/65028096 + (22231132110163*z^14)/520224768 - (86269923640667*z^13)/1040449536 + (659908638483491*z^12)/5549064192 - (614010908943583*z^11)/5549064192 + (1542961696333139*z^10)/66588770304 + (3684172462923767*z^9)/29595009024 - (86167575482678911*z^8)/355140108288 + (125425715847634333*z^7)/532710162432 - (428976223694222489*z^6)/4261681299456 - (115352974498514309*z^5)/2130840649728 + (2032084483563487*z^4)/16911433728 - (1527197305978447*z^3)/16911433728 + (778786651869991*z^2)/21743271936 - (134768480539*z)/18874368 + 2325457729/4194304, z, 1)^(1/2)' literally. Use 'isAlways' to test the statement mathematically.

Error in indexing (line 920)
X = find(mupadmex('symobj::logical',A.s,9)) - 1;

Error in indexing (line 968)
R_tilde = builtin('subsref',L_tilde,Idx);

disp(['q: ', num2str(length(positive_solutions))]);

Answer 1

[Dyuman Joshi]

solve() is unable to find an explicit solution to the given equation, which is expected given the degree of the polynomial - on simplifying the equation, you can observe that the highest power is 40.

There is no known method to find explicit solution for a polynomial of that degree.

You can use vpasolve -

Also, you should update the condition to check for non-complex roots as well.

syms r;

eq = (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^2 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^4 * sqrt(((1/2) + 1)^2 - r^2)^2 * r^4 ...
    + (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^8 ...
    - 2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) + 1)^2 - r^2)^2 * sqrt(((1/2) - 1)^2 - r^2)^2 * r^8 ...
    + (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 * r^8;

disp(simplify(eq))

sol = vpasolve(eq, r)

sol = 

%Update the condition
positive_solutions = sol(real(sol) > 0 & imag(sol) == 0 );
disp(['q: ', num2str(length(positive_solutions))]);

q: 2

[Imane Zemmouri]

Thank you very much.

[Dyuman Joshi]

You're welcome!

Answer 2

[Walter Roberson]

Your equality is equivalent to a polynomial of degree 20.

MATLAB is not able to find a closed form solution for the roots (which is not surprising -- theory says that degree 4 is the largest degree that can routinely be solved algebraically.) So MATLAB returns a form that "stands in" for the polynomial roots.

You then ask to reduce the solutions to the positive numbers. That requires being able to prove that particular solutions are positive... which is difficult when you cannot express the solutions explicitly.

[Walter Roberson]

syms r;

eq = (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^2 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^4 ...
    - 2 * (3040*r^6 - 12432*r^4 + 18544*r^2 - 18592)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) - 1)^2 - r^2)^4 * sqrt(((1/2) + 1)^2 - r^2)^2 * r^4 ...
    + (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^4 * sqrt(((1/2) + 1)^2 - r^2)^4 * r^8 ...
    - 2 * (768*r^8 - 768*r^6 + 288*r^4 - 304*r^2 + 69)^2 * (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^2 ...
    * sqrt(((1/2) + 1)^2 - r^2)^2 * sqrt(((1/2) - 1)^2 - r^2)^2 * r^8 ...
    + (768*r^8 - 6912*r^6 + 23328*r^4 - 35248*r^2 + 33081)^4 * sqrt(((1/2) - 1)^2 - r^2)^4 * r^8;

sol = solve(eq, r);
dsol = double(sol);
positive_solutions = sol(imag(dsol)==0 & real(dsol) > 0)

positive_solutions = 

disp(['q: ', num2str(length(positive_solutions))]);

q: 2

disp(double(positive_solutions))

    0.5000
    0.5000

[Imane Zemmouri]

Thank you very much.