Hello,
there is a precision problem in the numerical solutions found by algsys (or
solve). Here is an example run on maxima 5.17.1:
(%i1) algepsilon:10^16$
(%i2) fpprec:50$
(%i3) algsys([x^2+2*x*y^2,x*y+2*y^3-1],[y,x]),bfloat;
(%o3) [[y = 7.9370049732915826101731227026903070509433746337891b-1,
x = 2.7975867176063777611950563978988538844221110597592b-16],
[y = 6.8736481849930120002767353071249090135097503662109b-1 %i
- 3.9685026299204984034929566405480727553367614746094b-1,
x = - 4.1461656540724973579193248889167889387772229282426b-16 %i
- 4.0089773009749017276405892814572599004491524178964b-16],
[y = - 6.8736481849930120002767353071249090135097503662109b-1 %i
- 3.9685026299204984034929566405480727553367614746094b-1,
x = 4.1461656540724973579193248889167889387772229282426b-16 %i
- 4.0089773009749017276405892814572599004491524178964b-16]]
(%i4) algsys([x^2+2*x*y^2,x*y+2*y^3-1],[y,x]), algexact=true;
sqrt(3) %i - 1 sqrt(3) %i + 1
(%o4) [[y = --------------, x = 0], [y = - --------------, x = 0],
1/3 1/3
2 2 2 2
1
[y = ----, x =
0]]
1/3
2
(%i5) %,bfloat;
(%o5) [[y = 3.9685026299204986868792640981807706509787333197496b-1
(1.7320508075688772935274463415058723669428052538104b0 %i - 1.0b0),
x = 0.0b0], [y = - 3.9685026299204986868792640981807706509787333197496b-1
(1.7320508075688772935274463415058723669428052538104b0 %i + 1.0b0),
x = 0.0b0], [y = 7.9370052598409973737585281963615413019574666394993b-1,
x = 0.0b0]]
As you can see,even when increasing the default algepsilon and using big
floats, the real root y becomes incorrect very fast, since the correct
value is y = 7.937005259... while algsys (or solve) get y = 7.93700497..
hence there are only 6 or 7 correct figures. This is on a problem without
apparent numerical instabilities, without very large or very small
parameters, etc. in which case one may fear totally incoherent results.
This looks like a bug in the numerical algorithm of algsys.
By the way, i came to this example through Grobner basis: here is a way to
find the simple exact solution.
(%i6) load(grobner);
Loading maxima-grobner $Revision: 1.5 $ $Date: 2008/05/05 08:47:28 $
(%o6) /usr/local/share/maxima/5.17.1/share/contrib/Grobner/grobner.lisp
(%i7) poly_reduced_grobner([x^2+2*x*y^2,x*y+2*y^3-1],[y,x]);
3
(%o7) [x, 2 y - 1]
Here it is apparent that common roots of the ideal are such that x = 0, y is
a cubic root of 1/2.
Since people were asking information on grobner basis, here is a web site
which describes sufficiently clearly what it is about, and how grobner
basis are computed:
http://www.geocities.com/famancin/
--
Michel Talon