From: dcons@world.std.com (Denis Constales) Subject: Re: Areas Date: Fri, 09 Apr 1999 10:50:58 +0200 Newsgroups: sci.math Keywords: Many ways to compute the area of a circle In article <370D9D70.BDD98790@ns.sympatico.ca>, sylvester wrote: > What is the area of a circle given 3 non-collinear points ? There are several ways of answering this question. One would be to compute the sides a,b,c of the triangle defined by the three points, to apply the formula for the radius of the circumcircle, a b c R = 1/4 ------------------------------- sqrt(s (s - a) (s - b) (s - c)) where s = (a + b + c)/2 and to square it and multiply by Pi. Since a = sqrt((x2 - x3)^2 + (y2 - y3)^2) etc., there's lots of square roots occurring, and it takes somes simplificatory squirreling to obtain finally a result that does not contain square roots anymore. Another solution uses only coordinates, considering the system (x1-x0)^2+(y1-y0)^2=R^2 (x2-x0)^2+(y2-y0)^2=R^2 (x3-x0)^2+(y3-y0)^2=R^2 Expanding the squares gives rise to a linear system in the unknowns x0, y0, and (R^2 - x0^2 - y0^2): x1^2 + y1^2 = 2 * x1 * x0 + 2 * y1 * y0 + (R^2 - x0^2 - y0^2) etc. which can be solved using Cramer's method. From this, R^2 can be reconstructed, and the area is then simply Pi times this result. Summarizing, it is useful to introduce the following intermediate values: let D0 be the determinant of [1 x1 y1] [ ] [1 x2 y2] [ ] [1 x3 y3] D1 the determinant of [ 2 2 ] [x1 + y1 x1 y1] [ ] [ 2 2 ] [x2 + y2 x2 y2] [ ] [ 2 2 ] [x3 + y3 x3 y3] Dx the determinant of [ 2 2 ] [1 x1 + y1 y1] [ ] [ 2 2 ] [1 x2 + y2 y2] [ ] [ 2 2 ] [1 x3 + y3 y3] and Dy the determinant of [ 2 2] [1 x1 x1 + y1 ] [ ] [ 2 2] [1 x2 x2 + y2 ] [ ] [ 2 2] [1 x3 x3 + y3 ] Then the area is given by 2 2 Pi (Dx + Dy + 4 D0 D1) 1/4 ------------------------ 2 D0 Since the D's are polynomials in the coordinates, this is a rational function of them. For completeness, these are Maple procs for both approaches: way1 := proc(x1, y1, x2, y2, x3, y3) local a, b, c, s, R; a := sqrt((x2 - x3)^2 + (y2 - y3)^2); b := sqrt((x1 - x3)^2 + (y1 - y3)^2); c := sqrt((x2 - x1)^2 + (y2 - y1)^2); s := 1/2*a + 1/2*b + 1/2*c; R := 1/4*a*b*c/sqrt(s*(s - a)*(s - b)*(s - c)); Pi*R^2 end; way2 := proc(x1, y1, x2, y2, x3, y3) local D0, D1, Dx, Dy; D0 := x1*y2 - x1*y3 - y1*x2 + y1*x3 - x3*y2 + y3*x2; D1 := x1^2*y3*x2 - x1^2*x3*y2 + y1^2*y3*x2 - y1^2*x3*y2 - x2^2*x1*y3 + x2^2*y1*x3 - y2^2*x1*y3 + y2^2*y1*x3 + x3^2*x1*y2 - x3^2*y1*x2 + y3^2*x1*y2 - y3^2*y1*x2; Dx := x2^2*y3 + y2^2*y3 - y2*x3^2 - y2*y3^2 - y3*x1^2 - y3*y1^2 + x3^2*y1 + y3^2*y1 + y2*x1^2 + y2*y1^2 - x2^2*y1 - y2^2*y1; Dy := - x1^2*x2 + x1^2*x3 + x1*x2^2 + x1*y2^2 - x1*x3^2 - x1*y3^2 - y1^2*x2 + y1^2*x3 + x3^2*x2 - x3*x2^2 - x3*y2^2 + y3^2*x2; 1/4*Pi*(Dx^2 + Dy^2 + 4*D0*D1)/D0^2 end; Finally, note that when using the second method, it is numerically advantageous to subtract the coordinates (x1+x2+x3)/3 and (y1+y2+y3)/3 from all x resp. y coordinates first, to avoid cancellation of large terms in the determinants. -- Denis Constales - dcons@world.std.com - http://cage.rug.ac.be/~dc/ ============================================================================== From: dcons@world.std.com (Denis Constales) Subject: Re: Areas Date: Fri, 09 Apr 1999 14:15:41 +0200 Newsgroups: sci.math In article , dcons@world.std.com (Denis Constales) wrote: > [Snip!] Or, combining both methods in a third one, the area is given by 2 2 2 2 2 2 Pi ((x2 - x3) + (y2 - y3) ) ((x1 - x3) + (y1 - y3) ) ((x2 - x1) +(y2 - y1) ) 1/4-------------------------------------------------------------------------------- 2 (x1 y2 - x1 y3 - y1 x2 + y1 x3 - x3 y2 + y3 x2) or equivalently 2 2 2 2 2 2 Pi ((x2 - x3) + (y2 - y3) ) ((x1 - x3) + (y1 - y3) ) ((x2 - x1) + (y2 - y1) ) 1/4 -------------------------------------------------------------------------------- 2 ((x1 - x3) (y2 - y3) - (y1 - y3) (x2 - x3)) (this uses only differences of coordinates, which is better from a numerical viewpoint). Program: way3 := proc(x1, y1, x2, y2, x3, y3) local a2, b2, c2, D0; a2 := (x2 - x3)^2 + (y2 - y3)^2; b2 := (x1 - x3)^2 + (y1 - y3)^2; c2 := (x2 - x1)^2 + (y2 - y1)^2; D0 := (x1 - x3)*(y2 - y3) - (y1 - y3)*(x2 - x3); 1/4*Pi*a2*b2*c2/D0^2 end -- Denis Constales - dcons@world.std.com - http://cage.rug.ac.be/~dc/