Singular

Thank you. I really appreciate your help. I'll be honest with you, before last week I didn't know anything (zero) about "rings" and "ideals". Also, I am not a mathematician, I am an experimental physicist.
I just have, on an event by event base, sets of polynomials with 3 unknowns and with real coefficients (no rationals) and found out that in order to solve them I need either homotopy or this whole ring stuff. In my case I know the polynomial sets always have at least one solution, because they are the result of a physics process with at least one intersection point.
I am currently reading the book "Ideals, Varieties, and Algorithms" by Cox et al. to learn the basics of all that. So, I think it's better for me to read this book and just use the code as you suggested to me (which seems to work very well) and hopefully the book will help me understanding why.
More or less, I am guessing I am first to solve an "approximate" set based on rationals, before going on to the real set using reals. But going on asking you about this wouldn't be really fair, so I will study the book until I can make sense of it.

Best regards,

Machiel

It is not true that three equations in three unknows always define a finite set of solutions.
It already fails for three planes in three dimensional space when the planes meet in a common axis.

To see the difference between the rings, use Singular as a calculator and do some simple experiments:
[code]
ring rr=real,x,dp;
2/3+1/5
ring r0 = 0,x,dp;
2/3+1/5;
[/code]
As I explained in my first answer, the problem is that the standard basis computation over the
reals is not stable. It "destroys" the defining equation. Call std(i_ideal); to see what I mean.

If you work over the rationals, then the coefficients are always exact fractions and there is no
problem with inaccuray. Numerical solving, i.e. compuations over the reals ,should be the last
not the first step.

-----
P.S.: With the[b] User Control Panel [/b] of this forum it should be possible to configure automatic notification.

I am not familiar with Maple.
However, in your example, you put
[code]ring r0 = 0,(x,y,z),dp;[/code]
According to the manual, this means a ring of "rationals".
Indeed I have to put this line, because without it, the code complains that the ideal is not zero-dimensional. (Which does not make sense either, because a set of 3 equations with 3 unknowns actually is zero-dimensional, right?).
What I don't see is why I have to create a ring of rationals, given that both the solutions as the coordinates of the polynomials are not rationals but reals.

By the way (unrelated), how can I configure this forum so that it sends me an email automatically when there is an answer in this thread?

Consult Singular's online manual, start with [b]2.3 Getting started [/b]
[url]http://www.singular.uni-kl.de/Manual/latest/sing_5.htm#SEC11[/url]
and follow the links given there.

Some hints o clearify what is going on there:
In computations with Maple, you have to specify at each command what the variables are,
but Singular uses a different concept.

You specify a ring, as a first idea it is a polynomial ring over a field,
by defining the characteristic of the coefficient, the names of the variables and the ordering of the terms.
You give the ring a name. e.g. r0, and this is then [b]basering[/b] you are working in.

With this definition you working so to speak in a "categeory". You can defined several
rings, switch between them by calling [b]setring[/b] and transfer the objects
which can be [b]ideals[/b], [b]numbers[/b], [b]matrices[/b] etc. by calling [b]imap[/b], [b]fetch[/b], [b]map[/b].

Hi gorzel,

Thank you very much, following your advice I now do this:

[code]
ring s = real,(x,y,z),dp;

poly f = 723.056*x^2 + 12470.8*x*y - 86766.7*x + 7923.06*y^2 - 748246*y - 2876.94*z^2 + 8630.83*z + 2.59653e+06;
poly g = 10806.8*x^2 + 648407*x - 3593.22*y^2 + 373419*y - 3593.22*z^2 + 10779.7*z + 16315;
poly h = -782.305*x^2 - 12470.8*x*y - 694938*x + 6417.69*y^2 + 292823*y - 4382.31*z^2 + 13146.9*z - 2.82616e+06;

ideal i_ideal=f,g,h;

ring r0 = 0,(x,y,z),dp;
ideal I = imap(s,i_ideal);

LIB "solve.lib";
solve(I);
[/code]

Which gives the expected results.
The only problem is that I don't understand what I am doing. I understand that a "ring" is a collection of polynomials, and since I work with real numbers I supposed I had to define a ring of real numbers.
But what exactly is "ring r0" doing and why is it defined as "0"? It's never explicitly called later anymore. And what is an "imap"? And what does it mean to create an "ideal" from another ideal and a map s? And why do you say that "my ideal is reduced to the constant 1"?
I am aware that these are probably very fundamental questions, but whenever I look for answers I either find very theoretical mathematical papers or Wikipedia but neither really can explain to me how to relate those terms to the basic problem of three polynomials of degree 3 with 3 unknowns, which is the problem I am interesting in. So, if you can forward me to some instructive page/article where these concepts ("ideal" "ring" etc..) are introduced for a beginner, I would be very grateful.
But in any case, with your help for now I am already very happy, at least now I can solve my polynomials.

Regards,

Machiel

Standard basis computations over the reals are not numerical stable.
[b]Solve[/b] makes internal a check [code] dim(std(i))==0[/code]
but in your example the ideal is reduced to the constant 1, defining the empty set with dimension -1.
At that point [b]solve[/b] stopped with an ERROR message.

I the second example you gave the input has rational (integer) coefficients
but in your post it is not indicated what kind of basering you have chosen.

It may work over the reals, but you get preciser (integer) results when you call it within
[code]ring r =0,(x,y,z),dp; [/code]
General rules how to use solve:
[list]
1. Work over the Q[x,y,z]

2. if (size(variables(I)) < nvars(basering)) then reduce the variables in the defined ring

3. if (dim(std(I)) > 0) then make a primary decomposition (or call facstd) to separate the
higher dimensional components.

4. If (dim(std(I)==0) then vdim(std(I)) gives the number of solutions counted with multplicities

5. call solve(I,[...]).
[/list]
You could fetch your input to the rationals
(but if the input were rounded numbers, the exact rationals are not precise)
[code]
> ring r0 = 0,(x,y,z),dp;
> ideal I = imap(s,i_ideal);
> dim(std(I));
0
> vdim(std(I));
8
> solve(I);
// eight solutions, two of them are complex
[/code]
A naive method to filter the complex solutions is to check whether
the parameter [i]i[/i] occurs:
[code]
for (int j=1;j<=ncols(SOL);j++) {
if (!find(string(SOL[j]),"i)) { SOL[j]; }
}
[/code]

Hello,

I execute the following code:

[code]ring s = real,(x,y,z),dp;
poly f = 723.056*x^2 + 12470.8*x*y - 86766.7*x + 7923.06*y^2 - 748246*y - 2876.94*z^2 + 8630.83*z + 2.59653e+06;
poly g = 10806.8*x^2 + 648407*x - 3593.22*y^2 + 373419*y - 3593.22*z^2 + 10779.7*z + 16315;
poly h = -782.305*x^2 - 12470.8*x*y - 694938*x + 6417.69*y^2 + 292823*y - 4382.31*z^2 + 13146.9*z - 2.82616e+06;
ideal i_ideal=f,g,h;
option(redSB);
LIB "solve.lib";
def T=solve(i_ideal, 5, 0, "nodisplay");
[/code]

I get the error message:
[quote]? ideal not zero-dimensional[/quote]

Why? What does that mean? If I try a different set of polynomials, e.g.
[quote]poly f = 4*x^2 + 2*x*y + 3*x*z + 2*y*z + y^2 + z^2 - 42;
poly g = 6*x^2 + x*y - x*z - y*z + y^2 + z^2 - 12;
poly h = x^2 - x*y - 2*y^2 - 3*z^2 + 36;
[/quote]
the script works fine. However, with the previous set and different software, I do find solutions, so it's not a set without solutions.

By the way, another question is: when I want to print the solutions, I could only find the method
[code]setring T;
SOL;
[/code]
However, this gives all complex solutions. Is there a way to print out only the real solutions?

Many thanks,

Machiel