2
$\begingroup$

I have two dynamic perpendicular parabolas: $y=ax^2+c, x=by^2+d$, which intersect in four points (the coefficients are chosen that way). There is a fact that we can construct a circle through these four points.

My question is which command allows to construct a circle through four points?

$\endgroup$
4
  • 1
    $\begingroup$ Three points determine a circle. $\endgroup$
    – Mark Adler
    Commented May 15, 2016 at 20:00
  • $\begingroup$ yes, and i can graph it. But because of the machine accuracy being nonezero it can pass near the fourth point $\endgroup$ Commented May 15, 2016 at 20:19
  • 1
    $\begingroup$ @Michael Freimann If you edit your post to include the smallest and nicest (indented four spaces to let SE know this is code and not to be desktop published) bit of code, NOT a screenshot image, that shows your two parabolas with sliders then someone else might take a few minutes to try to write a couple of lines of code that would try to fit the best circle onto your four points. $\endgroup$
    – Bill
    Commented May 15, 2016 at 20:48
  • 1
    $\begingroup$ "There is a fact that we can construct a circle through these four points." - there's a constructive proof of that statement; why not use Mathematica to work through it? $\endgroup$ Commented May 16, 2016 at 2:29

2 Answers 2

17
$\begingroup$

This is a motivating post. I have made no effort to deal with intersection <4 pts. I post only illustrative examples.

The centre of the circle is the intersection of perpendicular bisectors of chords.

f[a_, b_, c_, d_] := 
 Module[{e1 = a x^2 + b, e2 = c y^2 + d, s, p, l, ln, ctr},
  s = Solve[{x, e1} == {e2, y}, {x, y}, Reals];
  p = {#, a #^2 + b} & /@ N[x /. s];
  l = Partition[p, 2, 1];
  ln = InfiniteLine[{(#1 + #2)/2, (#1 + #2)/2 + Cross[#2 - #1]}] & @@@
     l;
  ctr = RegionIntersection @@ ln[[{1, 2}]];
  ParametricPlot[{{x, e1}, {e2, y}}, {x, -3, 3}, {y, -3, 3}, 
   Epilog -> {Red, PointSize[0.03], Point[p], Black, Dashed, ln, 
     Green, ctr, Purple, Thick, 
     Circle[ctr[[1]], Norm[ctr[[1]] - p[[1]]]]}, 
   AspectRatio -> Automatic, 
   PlotLabel -> Row[{"{", e1, ",", e2, "}"}]]]

Examples:

Grid[Partition[
  f @@@ {{1, -2, 1, -2}, {1, -3, -1, 3}, {1, -4, 1, -2}, {-1, 2, 
     1, -3}}, 2]]

enter image description here

I reiterate I have made no effort to deal with number of intersections <4.

$\endgroup$
3
  • $\begingroup$ ...and of course, Circumsphere[] is built-in, and works in 2D. :) $\endgroup$ Commented May 16, 2016 at 9:31
  • $\begingroup$ @J.M. Thanks...didn't know that...always learning. :) $\endgroup$
    – ubpdqn
    Commented May 16, 2016 at 9:34
  • $\begingroup$ It's a convenient function, but I think your approach is more pedagogically appropriate for the OP. $\endgroup$ Commented May 16, 2016 at 9:45
6
$\begingroup$

Here is a standard algebraic approach, to complement @ubpdqn's geometric one.

If we have a sequence of points $(x_i,y_i)$, $i=1,2,\dots,n$ lying on a circle, then the linear system $$A\,(x_i^2+y_i^2)+B\,x_i+C\,y_i+D=0,\quad i=1,2,\dots,n$$ will have a nontrivial solution $(A_0,B_0,C_0,D_0)$. If there are at least three distinct points not on a line, the solution space will have dimension 1 and the circle through the points will be unique.

Here is a function that will find the intersections of curves and see a circle runs through the points.

ClearAll[circleIntersect];
circleIntersect::nocirc = "There is not a unique circle through the points.";
circleIntersect[eqs_, vars : {_, _}] := 
 Module[{circle, mat, form, res},
  form = {vars.vars, Sequence @@ vars, 1};       (* {x^2+y^2, x, y, 1} *)
  mat = form /. Solve[eqs, vars, Reals];         (* NSolve is often much faster *)
  Check[                   (* Dot::dotsh, Set::shape errors if null space not 1D *)
   {circle} = Sow[NullSpace[mat], "basis"].form; (* make basis available via Reap *)
   res = circle == 0,
   Message[circleIntersect::nocirc];
   res = $Failed
   ];
  Sow[mat[[All, {2, 3}]], "points"];             (* make points available via Reap *)
  res /; FreeQ[res, $Failed]
  ]

Examples

Two parabolas intersecting in four points.

p1 = y == x^2 - 4;
p2 = x == 2 y^2 - 3;
{eq, {pts}} = Reap[circleIntersect[{p1, p2}, {x, y}], "points"];
If[FreeQ[eq, circleIntersect],
 ContourPlot[Evaluate@{p1, p2, eq}, {x, -4, 4}, {y, -4, 4},
  Epilog -> {Red, PointSize[Medium], Point@pts}]]

Mathematica graphics

Solve will return points in terms of Root objects. If there is a nice form for the solution, FullSimplify might find it:

eq // FullSimplify
(*  2 (x^2 + y^2) == 11 + x + 2 y  *)

Tangent parabolas require some care. In this case if the machine real 4. is replaced by the exact integer 4, Solve runs for longer than I'm willing to wait. Numeric error with the machine-precision 4. leads to nonzero imaginary parts in the solution; hence the use of Chop below. One could replace 4. by the arbitrary-precision 4.`16, which would lead to a zero imaginary part.

p1 = y == x^2 - 4.;
p2 = x == y^2 - Root[16411 - 1152 #1 - 4096 #1^2 + 256 #1^3 &, 2];
{eq, {pts}} = Reap[circleIntersect[{p1, p2}, {x, y}], "points"];
If[FreeQ[eq, circleIntersect],
 ContourPlot[Evaluate@{p1, p2, Chop@eq}, {x, -4, 4}, {y, -4, 4},
  Epilog -> {Red, PointSize[Medium], Point@pts}]]

Mathematica graphics

eq
(*
  (0.960573 + 9.29229*10^-16 I) + (0.16052 + 1.32742*10^-14 I) x +
     (0.16052 - 1.06369*10^-14 I) y - (0.16052 - 6.13072*10^-15 I) (x^2 + y^2) == 0
*)

Points not lying on a circle. This intersection of a parabola and a quartic does not define a circle.

p1 = y == x^2 - 4 - x^4/16;
p2 = x == 2 y^2 - 3;
circleIntersect[{p1, p2}, {x, y}]
If[FreeQ[%, circleIntersect],
 ContourPlot[Evaluate@{p1, p2, %}, {x, -4, 4}, {y, -4, 4}]]

Mathematica graphics

More than one circle. This intersection of two parabolas has only two points and does not define a unique circle.

p1 = y == x^2 - 4;
p2 = x == 2 y^2 - 1;
{eq, {{{pts}}, {{basis}}}} = 
  Reap[circleIntersect[{p1, p2}, {x, y}], {"points", "basis"}];
pts // N

Mathematica graphics

Note, one could detect this. The basis returned NullSpace consists of two vectors. Linear combinations of these parameterize all the circles through the two points of intersection.

Show[
 ContourPlot[Evaluate@{p1, p2}, {x, -4, 4}, {y, -4, 4}, ContourStyle -> Black],
 ContourPlot[  (* family of circles *)
  Evaluate[
   Table[{(1 - t^2)/(1 + t^2), (2 t)/(1 + t^2)}, {t, 1.8^Range[-8, 8]}] . 
    N@basis.(List /@ {x^2 + y^2, x, y, 1})], 
  {x, -4, 4}, {y, -4, 4},
  Contours -> {{0}}]
 ]

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ beautiful...+1 of course :) $\endgroup$
    – ubpdqn
    Commented May 17, 2016 at 2:47

Not the answer you're looking for? Browse other questions tagged or ask your own question.