Menu
▾
▴
#3013 'solve' can depend on order (or presence) of variables given to solve
Problem with the 'solve' function - results/whether or not it returns anything dependent depend on the order (or presense) of the list of variables given to solve -- consider following MWE:
like : 28*log((1-s)*r)+18*log(s*(1-s)*r)+18*log(s*s*(1-s)*r)
+436*log((1-s)*(1-r)+s*(1-s)*(1-r)+s*s*(1-s)*(1-r)+s*s*s);
dlike_s : diff(like,s)=0;
dlike_r : diff(like,r)=0;
If I try solve([dlike_s,dlike_r],[s,r]);
(%i4) solve([dlike_s,dlike_r],[s,r]);
(%o4) []
if I simply switch the order of the variables...
(%i5) solve([dlike_s,dlike_r],[r,s]);
(%o5) [[r = 50653/(454750+4000*sqrt(4021)),s = -(sqrt(4021)+5)/74],
[r = -50653/(4000*sqrt(4021)-454750),s = (sqrt(4021)-5)/74]]
Further, if I leave the list of variables off altogether:
(%i6) solve([dlike_s,dlike_r]);
(%o6) [[r = 50653/(454750+4000*sqrt(4021)),s = -(sqrt(4021)+5)/74],
[r = -50653/(4000*sqrt(4021)-454750),s = (sqrt(4021)-5)/74]]
Also, same is true if 'algsys' is used instead of 'solve'.
Diff:
Diff:
Another example from the mailing list 2016-05-11: "Solve command succeeds or fails based on order variables are specified in the argument"
The first is definitely an algsys problem. I will have a look.
I have simplified the original equations and changed variables, but the bug is still there.
I have traced the working solution through algsys and followed the process manually in maxima. Here is the manual workings for the first solution found by maxima at %o3 (see batch file bug3013.mac attached). It is the second solution in the list.
Last edit: David Billinghurst 2016-10-05
This is what maxima attempts for the failed case.
1. algsys choses to eliminate y first
2. we solve for x successfully
3. substitute x into eq1 and eq2
4. equations for y are cubic
5. solutions are correct, but too complicated as maxima can't denest cube roots
6. fail
We could make progress by any of:
1. chosing to first eliminate x. This is based on a heuristic and it may be possible to do better.
2. denesting cube roots. Not trivial but the correct solution.
3. switching to a numerical solution at some point. This is done if solve can't find a symbolic solution, but not when solve finds a "complicated" solution. Perhaps we could numerically evaluate constant expressions in presultant and do something - at least give a warning - if the result is close to zero.
(%i2) eq1:125*x*y^3-125*x+16 3 (%o2) 125 x y - 125 x + 16 (%i3) eq2:713*x*y^4-681*x*y^3+(-59)*x*y+59*y+27*x-27 4 3 (%o3) 713 x y - 681 x y - 59 x y + 59 y + 27 x - 27 /* Choose to eliminate y - (leastvars eq1 eq2 ) => eq1 y */ (%i4) r:resultant(eq1,eq2,y) 3 2 (%o4) 20720464 x (2812500 x - 909500 x + 50653) /* Factor resultant Discard conants and duplicates => rf1 and rf2 */ (%i5) rf:factor(r) 3 2 (%o5) 20720464 x (2812500 x - 909500 x + 50653) (%i6) rf1:inpart(rf,2) 3 (%o6) x (%i7) rf1:inpart(rf,1) (%o7) 20720464 (%i8) rf2:inpart(rf,3) 2 (%o8) 2812500 x - 909500 x + 50653 /* solve rf1 and rf2 for x X is ordered to match algsys */ (%i9) X1:solve(rf1,x) (%o9) [] (%i10) ev(X1:ratsimp(X1),algebraic = true) (%o10) [] (%i11) X2:solve(rf2,x) 16 sqrt(4021) - 1819 16 sqrt(4021) + 1819 (%o11) [x = - --------------------, x = --------------------] 11250 11250 (%i12) ev(X2:ratsimp(X2),algebraic = true) 16 sqrt(4021) - 1819 16 sqrt(4021) + 1819 (%o12) [x = - --------------------, x = --------------------] 11250 11250 (%i13) X:append(X1,reverse(X2)) 16 sqrt(4021) + 1819 16 sqrt(4021) - 1819 (%o13) [x = --------------------, x = - --------------------] 11250 11250 /* Subst X[1] into eq1 and eq2 Confirm it isn't a solution Don't reproduce. */ /* Now subst X[2] into eq1 and eq2 - solve for y */ (%i14) eq12:subst(X[2],eq1) 3 (16 sqrt(4021) - 1819) y 16 sqrt(4021) - 1819 (%o14) (- -------------------------) + -------------------- + 16 90 90 (%i15) eq22:subst(X[2],eq2) 4 3 713 (16 sqrt(4021) - 1819) y 227 (16 sqrt(4021) - 1819) y (%o15) (- -----------------------------) + ----------------------------- 11250 3750 59 (16 sqrt(4021) - 1819) y 3 (16 sqrt(4021) - 1819) + --------------------------- + 59 y - ------------------------ - 27 11250 1250 /* Choose to work on eq12 OK. degree(eq12)=3 < degree(eq22) = 4 */ (%i16) ev(eq12:ratsimp(eq12*90),algebraic = true) 3 (%o16) (1819 - 16 sqrt(4021)) y + 16 sqrt(4021) - 379 (%i17) ev(eq22:ratsimp(eq22*11250),algebraic = true) 4 3 (%o17) (1296947 - 11408 sqrt(4021)) y + (10896 sqrt(4021) - 1238739) y + (944 sqrt(4021) + 556429) y - 432 sqrt(4021) - 254637 (%i18) r2:resultant(eq22,eq12,y) (%o18) 0 /* = 0. Correct */ (%i19) Y2:solve(eq12,y) 1/3 1/3 sqrt(3) (16 sqrt(4021) - 379) %i - (16 sqrt(4021) - 379) (%o19) [y = --------------------------------------------------------------, 1/3 2 (16 sqrt(4021) - 1819) 1/3 1/3 sqrt(3) (16 sqrt(4021) - 379) %i + (16 sqrt(4021) - 379) y = - --------------------------------------------------------------, 1/3 2 (16 sqrt(4021) - 1819) 1/3 (16 sqrt(4021) - 379) y = -------------------------] 1/3 (16 sqrt(4021) - 1819) /* Here is the failure. The third root is solution we want, but we can't denest cube root Later we will reject this solution */ (%i20) rectform(float(Y2)) (%o20) [y = 0.4622388881753429 - 0.8006212395538426 %i, y = 0.8006212395538426 %i + 0.4622388881753429, y = - 0.9244777763506857] (%i21) float(-(sqrt(4021)+5)/74); (%o21) - 0.9244777763506861Last edit: David Billinghurst 2016-10-05
Here is a reduced test case for the failure above. It should return the real solution of Y2 above.