Re: Re: More about ellipse and circle intersection
- To: mathgroup at smc.vnet.net
- Subject: [mg38531] Re: [mg38439] Re: [mg38426] More about ellipse and circle intersection
- From: Tomas Garza <tgarza01 at prodigy.net.mx>
- Date: Sun, 22 Dec 2002 04:14:00 -0500 (EST)
- References: <200212140820.DAA09217@smc.vnet.net> <200212150710.CAA15940@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Definitely, Daniel's explanation solves the mistery (even though most of the
details go well above my head). Cf. his comment below:
> The reason is that the generic radical solution does not work for those
> particular values of parameters.
Those particular values of the parameters were chosen, if not randomly,
certainly with no particular intention in mind. They were meant only to
illustrate the method, and they just happened to create some kind of havoc.
The morale is that utmost care should be taken when working with tools such
as Solve. The main problem is that not only the generic solution does not
work always (such as in the example here), but that it may yield a wrong
numerical solution. It just happened that I could visually determine that
the solution was wrong because I could plot the situation. Otherwise I might
have been in deep trouble.
Tomas Garza
----- Original Message -----
From: "Daniel Lichtblau" <danl at wolfram.com>
To: mathgroup at smc.vnet.net
Subject: [mg38531] [mg38439] Re: [mg38426] More about ellipse and circle intersection
> Tomas Garza wrote:
> >
> > In a previous message [mg38381] Re: Ellipse and circle intersection, I
> > hinted that Solve would not always give the right answer when trying to
> > obtain the intersections of an ellipse and a circle. Two comments were
> > received (Tom Burton and Rasmus Debitsch) suggesting that there was
> > nothing wrong. Still, definitely there is, as shown in the example
> > below. The situation is as follows:
> >
> > In[1]:=
> > ellipse = (x - c)^2/b^2 + (y - d)^2/a^2 == 1;
> > circ = x^2 + y^2 == 1;
> > sol = Solve[{ellipse, circ}, {x, y}];
> >
> > The solution comes out all right in terms of the four parameters a, b,
> > c, d. The following sets of values for the parametrs are tested:
> >
> > In[2]:=
> > example1Tom = {a -> 1.2, b -> 1.3, c -> 0.2, d -> 0.3};
> > example2Tom = {a -> 1.2, b -> 1.3, c -> 1.2, d -> 1.3};
> > example3Rasmus = {a -> 2, b -> 1, c -> 1, d -> 1};
> > example4Tomas = {a -> 0.25, b -> 0.75, c -> 0.5, d -> 0};
> >
> > (the first two come from Tom Burton, the third one from Rasmus, and the
> > fourth one is mine). Numerical solutions are then obtained for each set:
> >
> > In[3]:=
> > sol1=sol/.example1Tom;
> > In[4]:=
> > sol2=sol/.example2Tom;
> > In[5]:=
> > sol3=sol/.example3Rasmus//N;
> > In[6]:=
> > sol4=sol/.example4Tomas;
> >
> > In each of the first three cases the correct intersections (as observed
> > graphically through ImplicitPlot) are found, in addition to some complex
> > points. However, the fourth case fails to give a correct answer, even if
> > the two curves intersect very nicely at four different points in the
> > plane. This points to a weird behavior of Solve. I hope someone will
> > come forward with some explanation.
> >
> > Tomas Garza
> > Mexico City
>
>
> The reason is that the generic radical solution does not work for those
> particular values of parameters. We can see why as follows. First we
> create polynomials and find a quartic in one variable.
>
> eqns = {ellipse,circ};
> polys = Map[#[[1]]-#[[2]]&, eqns];
>
> InputForm[uniploy = First[GroebnerBasis[polys, {x,y},
> CoefficientDomain->RationalFunctions]]]
>
> Out[67]//InputForm=
> 1 + b^(-4) - 2/b^2 - (2*c^2)/b^4 - (2*c^2)/b^2 + c^4/b^4 - (2*d^2)/a^2 +
> (2*d^2)/(a^2*b^2) + (2*c^2*d^2)/(a^2*b^2) + d^4/a^4 +
> ((4*d)/a^2 - (4*d)/(a^2*b^2) - (4*c^2*d)/(a^2*b^2) - (4*d^3)/a^4)*y +
> (-2/a^2 - 2/b^4 + 2/b^2 + 2/(a^2*b^2) + (2*c^2)/b^4 + (2*c^2)/(a^2*b^2)
> +
> (6*d^2)/a^4 - (2*d^2)/(a^2*b^2))*y^2 + ((-4*d)/a^4 +
> (4*d)/(a^2*b^2))*
> y^3 + (a^(-4) + b^(-4) - 2/(a^2*b^2))*y^4
>
> We have a quartic univariate in y (in terms of the problem parameters).
>
> Now we plug in your values:
>
> InputForm[unipoly /. {a->1/4, b->3/4, c->1/2, d->0}]
> Out[69]//InputForm= (81*(-5/3 + (1024*y^2)/27 + (16384*y^4)/81))/16384
>
> Notice that it is biquadratic (that is, quadratic in y^2). But this is
> exactly the case for which the Cardano formulation of the quartic
> radical solution does not work.
>
> To obtain a solution that does not have this limitation one must avoid a
> radical formulation. This can be done as below. The solution (not shown)
> is in terms of parametrized Root[] functions.
>
> SetOptions[Roots, Quartics->False];
> soln = Solve[{ellipse,circ}, {x,y}];
>
> In[73]:= InputForm[soln /. N[{a->1/4, b->3/4, c->1/2, d->0}]]
> Out[73]//InputForm=
> {{x -> 0.981455818030629, y -> -0.19168849014437103},
> {x -> 0.981455818030629, y -> 0.19168849014437103},
> {x -> -1.106455818030629 - 8.619810689307775*^-27*I,
> y -> 1.1376714700528284*^-27 - 0.473544588453747*I},
> {x -> -1.106455818030629 + 8.619810689307775*^-27*I,
> y -> 1.1376714700528284*^-27 + 0.473544588453747*I}}
>
> If you are curious as to why the Cardano method fails when the quartic
> reduces to a biquadratic, I can send a derivation of that formula.
> Actually I should have put it in
>
>
http://library.wolfram.com/conferences/conference98/abstracts/various_ways_t
o_tackle_algebraic_equations.html
>
> In the notebook accessible from that link I show how one might derive
> the cubic radical solution. My updated version of the notebook also
> covers the quartic derivation so I will try to get that version put at
> the web site. The essence of the nongeneric problem is that one is
> confronted with a hidden division by zero.
>
>
> Daniel Lichtblau
> Wolfram Research
>
>
- References:
- More about ellipse and circle intersection
- From: Tomas Garza <tgarza01@prodigy.net.mx>
- Re: More about ellipse and circle intersection
- From: Daniel Lichtblau <danl@wolfram.com>
- More about ellipse and circle intersection