Finding the centroid of the area between two curves
Edit: Preface
After being on this site for a few months, I've learned quite a bit about what Mathematica can do quickly. Below I've revised my code to take advantage of built-in vectorized functions. One could compile the functions as well, but the advantages perhaps outweigh the disadvantages. The increase in speed is slight on approximate real numbers (because the number crunching is done in a very few function calls and so is already "mostly compiled" anyway); however, the first function centroid1 is exact on polygons, and a compiled function would not return an exact result on exact input. I added a new one that implements a Simpson's Rule type integration that is exact on figures whose boundary can be parametrized by quadratic polynomials. I also updated the timings, since I have a new machine.
Revised original method
Originally I had done this for someone who had wanted to find centroids of polygons, for which NIntegrate seemed too slow. I thought it was pretty cool at the time and wanted to share it. I guess it's pretty clear that the simple answer to the original question of @a3f is that there is no built-in centroid function.
The code below replaces the previous code but is based on the same mathematical ideas as before. The line integrals around the boundary,
$$A={1\over2} \int x\;dy - y\;dx, \quad M_x=-\int x\,y\;dx, \quad M_y=\int x\,y\;dy\,,$
which by Green's Theorem represent the area, x moment, and y moment respectively, are approximated using a linear interpolation of the boundary at sample points.
area1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
Total[x RotateLeft[y] - RotateLeft[x] y]] / 2;
xMoment1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
Total[(x - RotateLeft[x]) (2 x y + y RotateLeft[x] +
x RotateLeft[y] + 2 RotateLeft[x] RotateLeft[y])]] / 6;
yMoment1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
Total[(RotateLeft[y] - y) (2 x y + y RotateLeft[x] +
x RotateLeft[y] + 2 RotateLeft[x] RotateLeft[y])]] / 6;
centroid1[poly_Polygon] := centroid1[First[poly]];
centroid1[poly_List?(Depth[#] == 4 &)] := (* multiple-polygon figure (e.g. Italy) *)
{#[[1]], #[[2]]} / #[[3]] &[
# . Sign@Last[#] & @ {xMoment1 /@ #, yMoment1 /@ #, area1 /@ #} & @ poly];
centroid1[poly_List?(Depth[#] == 3 &)] := {xMoment1[#], yMoment1[#]} / area1[#] & @ poly;
Polygons
As a test I loaded Italy and found its centroid:
itPoly = CountryData["Italy", "Polygon"];
centroid1[itPoly]
(* {12.0822, 42.7896} *)
It's the centroid of a "flat-earth" Italy, since no adjustment for spherical coordinates is made.
Graphics[{itPoly, PointSize[Large], Red, Point[centroid1[itPoly]]}]
The new code is faster:
centroid1[itPoly] // AbsoluteTiming (* new *)
(* {0.000468, {12.0822, 42.7896}} *)
oldCentroid[itPoly] // AbsoluteTiming (* old *)
(* {0.007083, {12.0822, 42.7896}} *)
The function centroid1 is faster than NIntegrate on polygons:
{#[[1]], #[[2]]}/#[[3]] &@(Plus @@ (Module[{X, Y, a, xm, ym},
X = Interpolation[Append[#, First[#]] &[First /@ #], InterpolationOrder -> 1];
Y = Interpolation[Append[#, First[#]] &[Last /@ #], InterpolationOrder -> 1];
a = NIntegrate[X[t] Y'[t], {t, 1, 1 + Length[#]}];
xm = NIntegrate[-X[t] Y[t] X'[t], {t, 1, 1 + Length[#]}];
ym = NIntegrate[X[t] Y[t] Y'[t], {t, 1, 1 + Length[#]}];
{xm, ym, a}
] & /@ First[itPoly])) // AbsoluteTiming
(* {2.616189, {12.0822, 42.7896}} *)
Curves
The function centroid1 above is not exact on curves, but it may be accurate enough for some work. In most cases, NIntegrate is hard to beat it for accuracy. To apply centroid1, first you have to convert the curve to an approximate polygon (pick sample points).
Graphs
For the example proposed by @a3f, we can extract the plotted points of the graph with Cases. Since the base is straight, they form the polygon and we can apply centroid1:
(paraPoly = Transpose@{#, 1 - #^2/4} &@Range[-2., 2, 4/1000];
centroid1[paraPoly]) // AbsoluteTiming
(* {0.000356, {6.83717*10^-17, 0.4}} *)
The error in the y-coordinate is less than $10^{-6}$.
Graphics[{Polygon[paraPoly], PointSize[Large], Red,
Point[centroid1[paraPoly]]}, Frame -> True]
If we compare with Integrate and NIntegrate as @murray, @george2079 have answered, we get:
(a = Integrate[(1 - x^2/4), {x, -2, 2}];
xm = Integrate[x (1 - x^2/4), {x, -2, 2}];
ym = Integrate[(1 - x^2/4)/2 (1 - x^2/4), {x, -2, 2}];
{xm, ym} / a) // AbsoluteTiming
(* {0.020392, {0, 2/5}} *)
(a = NIntegrate[(1 - x^2/4), {x, -2, 2}];
xm = NIntegrate[x (1 - x^2/4), {x, -2, 2}];
ym = NIntegrate[(1 - x^2/4)/2 (1 - x^2/4), {x, -2, 2}];
{xm, ym} / a) // AbsoluteTiming
(* {0.013740, {0., 0.4}} *)
The first is exact and the second is accurate to $MachinePrecision.
Parametric Curves
I made up, at random, simple closed curve (picture down below),
f[t_] := {Cos[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t]), Sin[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t])};
and found that for 200000 points, the centroid1 calculation took about as long as the NIntegrate. So we can compare the results.
(polarPoly = Transpose @ f[Range[0., 1, 1/200000]];
centroid1[polarPoly]) // AbsoluteTiming
(* {0.104392, {0.95, -1.38784*10^-14}} *)
Module[{a, xm, ym},
a = NIntegrate[Evaluate[f[t][[1]] f'[t][[2]]], {t, 0, 1}];
xm = NIntegrate[Evaluate[-f[t][[1]] f[t][[2]] f'[t][[1]]], {t, 0, 1}];
ym = NIntegrate[Evaluate[f[t][[1]] f[t][[2]] f'[t][[2]]], {t, 0, 1}];
{xm, ym}/a
] // AbsoluteTiming
(* {0.108268, {0.95, -1.11203*10^-16}} *)
Module[{a, xm, ym},
a = Integrate[Evaluate[f[t][[1]] f'[t][[2]]], {t, 0, 1}];
xm = Integrate[-Evaluate[f[t][[1]] f[t][[2]] f'[t][[1]]], {t, 0, 1}];
ym = Integrate[ Evaluate[f[t][[1]] f[t][[2]] f'[t][[2]]], {t, 0, 1}];
{xm, ym}/a
] // AbsoluteTiming
(* {7.968150, {19/20, 0}} *)
The errors in the x coordinate are a little over 10^-10 and 10^-16 (or $MachinePrecision) for centroid1 and NIntegrate respectively. (But on 5000 points centroid1 is 100 times faster and has an error in x of a little over 10^-7, which is good enough for plotting.)
Addition -- Simpson's Rule
If we have an even number of points, we can use a method similar to Simpson's Rule. Think of the points as alternating between "endpoints" and "midpoints," and then use quadratic interpolation to approximate the boundary of the region. This tends to be better on curved boundaries and worse on straight boundaries such as polygons. With polygons, we will get an inaccurate result unless we bisect the segments; so from the points of view of both accuracy and speed, it loses over the linear method of centroid1.
With[{
areaFn = Total[(-3 #1 - 4 #2 + #3) #4 + 4 (#1 - #3) #5 - (#1 - 4 #2 - 3 #3) #6]/6 &,
momentFn = Total[(-10 #1 - 6 #2 + #3) #4^2 + (6 #1 - 8 #2 + 2 #3) (#4 #5) +
8 (#1 - #3) #5^2 + (-#1 + #3) (#4 #6) + (-2 #1 + 8 #2 - 6 #3) (#5 #6) +
(-#1 + 6 #2 + 10 #3) #6^2]/30 &},
area2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
areaFn[x[[1 ;; -1 ;; 2]], x[[2 ;; -1 ;; 2]],
RotateLeft[x[[1 ;; -1 ;; 2]]], y[[1 ;; -1 ;; 2]],
y[[2 ;; -1 ;; 2]], RotateLeft[y[[1 ;; -1 ;; 2]]]]
];
yMoment2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
momentFn[x[[1 ;; -1 ;; 2]], x[[2 ;; -1 ;; 2]],
RotateLeft[x[[1 ;; -1 ;; 2]]], y[[1 ;; -1 ;; 2]],
y[[2 ;; -1 ;; 2]], RotateLeft[y[[1 ;; -1 ;; 2]]]]
];
xMoment2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
-momentFn[y[[1 ;; -1 ;; 2]], y[[2 ;; -1 ;; 2]],
RotateLeft[y[[1 ;; -1 ;; 2]]], x[[1 ;; -1 ;; 2]],
x[[2 ;; -1 ;; 2]], RotateLeft[x[[1 ;; -1 ;; 2]]]]
]
];
centroid2[poly_Polygon] := centroid2[First[poly]];
centroid2[poly_List?(Depth[#] == 4 &)] /; And @@ (EvenQ /@ Length /@ poly) :=
{#[[1]], #[[2]]} / #[[3]] &[#.(Sign@Last[#]) &[
{xMoment2 /@ #, yMoment2 /@ #, area2 /@ #} &@ poly]];
centroid2[poly_List?(Depth[#] == 3 &)] /; EvenQ @ Length[poly] :=
{xMoment2[#], yMoment2[#]} / area2[#] &@ poly;
It's exact on regions with boundaries that have quadratic parametrizations. Here care is taken to have the bottom segment bisected by the origin as midpoint in order to get the exact answer.
(paraPoly = Transpose@{#, 1 - #^2/4} &@ Range[-2, 2, 4/4] ~Join~ {{0, 0}};
centroid2[paraPoly]) // AbsoluteTiming
(* {0.000289, {0, 2/5}} *)
On the parametric curve, with 1000 sample points, the error in x is a little over 10^-7, the same as centroid1 on 5000 points.
f[t_] := {Cos[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t]),
Sin[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t])};
(polarPoly = Transpose @ f[Range[0., 1, 1/999]];
centroid2[polarPoly]) // AbsoluteTiming
(* {0.000714, {0.95, 5.44169*10^-10}} *)
Now in V10 we have ImplicitRegion and RegionCentroid to do this easily:
reg = ImplicitRegion[y > 0 && y <= 1 - x^2/4, {x, y}];
Then:
RegionCentroid[reg]
{0, 2/5}
This is inelegant for the specific example, but may be useful for more general cases:
boolf[x_, y_] := (y > 0 && y < 1 - x^2/4);
area = NIntegrate[If[ boolf[x, y] , 1, 0], {x, -20, 20}, {y, -20, 20}]
-> 2.66667
cx = NIntegrate[If[ boolf[x, y] , x, 0], {x, -20, 20}, {y, -20, 20}]/ area
-> ~10^-16
cy = NIntegrate[If[ boolf[x, y] , y, 0], {x, -20, 20}, {y, -20, 20}]/area
-> 0.4
The exact results for the example are by the way:
area = Integrate[( 1 - x^2/4 ) , {x, -2, 2}] -> 8/3
cy = Integrate[( 1 - x^2/4 )/2 ( 1 - x^2/4 ) , {x, -2, 2}] / area -> 2/5
cx = Integrate[ x ( 1 - x^2/4 ) , {x, -2, 2}] / area -> 0
Edit:
Further explanation: we construct a logical function of x,y which is True inside your region and False otherwise, then perform area integrals over the entire plane using the definitions of area and centroid, and relying on the integrands to be zero outside the region of interest.
As noted by Daniel in a comment that seems to have disappeared(?) you can use the analytic Integrate function over the entire plane (±∞):
fa[x_,y_]:= If[(y > 0 && y < 1 - x^2/4), 1, 0];
area = Integrate[fa[x, y], {x, -∞, ∞}, {y, -∞, ∞}]
-> 8/3
cx = Integrate[ x fa[x,y], {x, -∞, ∞}, {y, -∞, ∞}]/ area
-> 0
cy = Integrate[ y fa[x,y], {x, -∞, ∞}, {y, -∞, ∞}]/area
-> 2/5
This does need a bit of warning, Integrate[] fails for even slightly more complicated expressions, so NIntegrate[] is a bit safer (though it yields a numerical approximation).