Green’s theorem and section properties

In the last post, I presented a simple Python module with functions for calculating section properties of polygons. Now we’ll go through the derivations of the formulas used in those functions.

The basis for all the formulas is Green’s theorem, which is usually presented something like this:

[\oint\limits_C P\; dx + Q\; dy = \iint\limits_A \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right)\; dx dy]

where [P] and [Q] are functions of [x] and [y], [A] is the region over which the right integral is being evaluated, and [C] is the boundary of that region. The integral on the right is evaluated in accordance with the right-hand rule, i.e., counterclockwise for the usual orientation of the [x] and [y] axes.

Region and nomenclature

The section properties of interest are all area integrals. We’ll use Green’s theorem to turn them into boundary integrals and then evaluate those integrals using the coordinates of the polygon’s vertices.


This is the easiest one, but instead of going through the full derivation here, I’ll refer you to this excellent StackExchange page by apnorton and just hit the highlights.

  1. The area is defined

    [A = \iint\limits_Adx dy]

    and we’ll choose [P = 0] and [Q = x] as our Greens’ theorem functions. This gives us

    [A = \iint\limits_Adx dy = \oint\limits_C x\; dy]
  2. We break the polygonal boundary into a series of straight-line segments, each of which can be parameterized this way:

    [x = x_i + (x_{i+1} - x_i)\;t] [y = y_i + (y_{i+1} - y_i)\;t] [dy = (y_{i+1} - y_1)\; dt]

    where the [(x_i, y_i)] are the coordinates of the vertices.

  3. Plugging these equations into the integral, we get

    [A = \frac{1}{2} \sum_{i=0}^{n-1}\; (x_{i+1} + x_i)(y_{i+1} - y_i)]

A note on the indexing: The polygon has [n] vertices, which we’ll number from 0 to [n-1]. The last segment of the boundary goes from [(x_{n-1}, y_{n-1})] to [(x_0, y_0)]. To make this work with the equation, we’ll define [(x_n, y_n) = (x_0, y_0)].

Let’s compare this with the area function in the module:

def area(pts):
  'Area of cross-section.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  s = 0
  for i in range(len(pts) - 1):
    s += x[i]*y[i+1] - x[i+1]*y[i]
  return s/2

We start by checking the pts list to see if the starting and ending items match. If they don’t, we copy the starting item to the end to fit the indexing convention discussed above. We then initialize some variables and execute a loop, summing terms along the way. Rewriting the loop in mathematical terms, we get

[A = \frac{1}{2} \sum_{i=0}^{n-1}\; x_i y_{i+1} - x_{i+1} y_i]

This doesn’t look like the equation derived from Green’s theorem, does it? But it’s not too hard to see that they are equivalent. Expanding out the binomial product in the earlier equation gives

[x_{i+1} y_{i+1} - x_{i+1} y_i + x_i y_{i+1} - x_i y_i]

As we loop through all the values of [i] from 0 to [n-1], the leading term of one trip through the loop will cancel the trailing term of the next trip through the loop. Here’s an example for a triangle:

[\quad (x_1 y_1 - x_1 y_0 + x_0 y_1 - x_0 y_0 )] [+\; (x_2 y_2 - x_2 y_1 + x_1 y_2 - x_1 y_1 )] [+\; (x_0 y_0 - x_0 y_2 + x_2 y_0 - x_2 y_2 )]

After the cancellations, all that’s left are the inner terms, and that’s the formula used in the area function.

The cancellation doesn’t do much for us here, changing from two additions and one multiplication per loop to two multiplications and one addition per loop. But we’ll see this same sort of cancellation in the other section properties, and it will provide greater simplification in those.


The centroid is essentially the average position of the area. If a sheet of material of uniform thickness and density were cut into a shape, the centroid would be the center of gravity, the balance point, of that shape. The coordinates of the centroid are defined this way:

[x_c = \frac{1}{A} \iint\limits_A x\; dx dy] [y_c = \frac{1}{A} \iint\limits_A y\; dx dy]

Let’s derive the formula for [x_c] for a polygon; the derivation of the formula for [y_c] will be similar.

In applying Green’s theorem, we’ll take [P = 0] and [Q = \frac{1}{2} x^2]. Therefore,

[x_c = \frac{1}{A} \iint\limits_A x\; dx dy = \oint\limits_C \frac{1}{2} x^2\; dy]

Breaking the polygonal boundary into straight-line segments and using the same parametric equations as before, we get an integral that looks like this

[\int_0^1 \frac{1}{2} \left[ (x_i^2 + 2 x_i (x_{i+1} - x_i)\; t + (x_{i+1} - x_i)^2\; t^2\right] (y_{i+1} - y_i)\; dt]

for each segment. This integral evaluates to

[\frac{1}{6} (x_{i+1}^2 + x_{i+1} x_i + x_i^2)\; ( y_{i+1} - y_i)]

so our formula for the centroid is

[x_c = \frac{1}{6A} \sum_{i=0}^{n-1} (x_{i+1}^2 + x_{i+1} x_i + x_i^2)\; ( y_{i+1} - y_i)]

As we found in the formula for area, the leading and trailing terms in the expansion of this product cancel out as we loop through the sum, leaving us with

[x_c = \frac{1}{6A} \sum_{i=0}^{n-1} -y_i x_{i+1}^2 + y_{i+1}x_{i+1} x_i - y_i x_{i+1} x_i + y_{i+1} x_i^2]

This looks like a mess, but it can be factored into a more compact form:

[x_c = \frac{1}{6A} \sum_{i=0}^{n-1} (x_{i+1} + x_i)\;(x_i y_{i+1} - x_{i+1} y_i)]

The expression for the other centroidal coordinate is as you’d expect:

[y_c = \frac{1}{6A} \sum_{i=0}^{n-1} (y_{i+1} + y_i)\;(x_i y_{i+1} - x_{i+1} y_i)]

These are the formulas used in the centroid function.

def centroid(pts):
  'Location of centroid.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  sx = sy = 0
  a = area(pts)
  for i in range(len(pts) - 1):
    sx += (x[i] + x[i+1])*(x[i]*y[i+1] - x[i+1]*y[i])
    sy += (y[i] + y[i+1])*(x[i]*y[i+1] - x[i+1]*y[i])
  return sx/(6*a), sy/(6*a)

Moments and product of inertia

You may be familiar with moments and products of inertia from dynamics, where the terms are related to the distribution of mass in a body. The moments and product of inertia we’ll be talking about here—more properly called the second moments of area—are mathematically similar and refer to the distribution of area across a planar shape.

The moments and product of inertia that matter in beam bending are taken about the centroidal axis (i.e., a set of [x] and [y] axes with the origin at the centroid of the shape). Since we don’t know where the centroid is when we set up our coordinate system, our list of vertex points don’t work off that basis. But we can still calculate the centroidal moments and product of inertia by using these formulas:

[I_{xx} = \iint\limits_A (y - y_c)^2\; dx dy] [I_{yy} = \iint\limits _A (x - x_c)^2 \; dx dy] [I_{xy} = \iint\limits_A (y - y_c)\; (x - x_c)\; dx dy]

We’ll concentrate on [I_{yy}]; the other two will be similarly derived.

First, let’s expand the square inside the integral and see what we get:

[I_{yy} = \iint\limits_A x^2\; dx dy - 2x_c \iint\limits_A x\; dx dy + x_c^2 \iint\limits_A dx dy]

The integral in the second term is [A x_c] and the integral in the third term is just [A]. Putting this together, we get1

[I_{yy} = \iint\limits_A x^2 dx dy - A\;x_c^2]

Since we already have formulas for [x] and [x_c], we can concentrate on the integral in the first term on the right.

Returning to Green’s theorem, we’ll use [P = 0] and [Q = \frac{1}{3}x^3], giving us

[\iint\limits_A x^2 dx dy = \oint\limits_C \frac{1}{3}x^3 \; dy]

Once again, we break the polygonal boundary into straight-line segments and use parametric equations to define the segments. For each segment, we’ll get the following integral:

[\int_0^1 \frac{1}{3} \left[ x_i^3 + 3 x_i^2 (x_{i+1} - x_1)\; t + 3 x_i (x_{i+1} - x_1)^2\; t^2 + (x_{i+1} - x_1)^3\; t^3 \right] (y_{i+1} - y_i) \; dt]

This integral evaluates to

[\frac{1}{12} \left[ x_{i+1}^3 + x_{i+1}^2 x_i + x_{i+1} x_i^2 + x_i^3 \right] (y_{i+1} - y_i)]

giving us

[\iint\limits_A x^2 dx dy = \frac{1}{12} \sum_{i=0}^{n-1} \left[ x_{i+1}^3 + x_{i+1}^2 x_i + x_{i+1} x_i^2 + x_i^3 \right] (y_{i+1} - y_i)]

Once again, if we expand out the product inside the sum, we’ll find that the leading and trailing terms cancel as we work through the loop. That gives us

[\frac{1}{12} \sum_{i=0}^{n-1} -y_i x_{i+1}^3 + y_{i+1} x_{i+1}^2 x_i - y_i x_{i+1}^2 x_i + y_{i+1} x_{i+1} x_i^2 - y_i x_{i+1} x_i^2 + y_{i+1} x_i^3]

And that long expression can be factored, leaving

[\iint\limits_A x^2 dx dy = \frac{1}{12} \sum_{i=0}^{n-1} (x_{i+1}^2 + x_{i+1} x_i + x_i^2)\; (x_i y_{i+1} - x_{i+1} y_i)]

Similar2 derivations give us

[\iint\limits_A y^2 dx dy = \frac{1}{12} \sum_{i=0}^{n-1} (y_{i+1}^2 + y_{i+1} y_i + y_i^2)\; (x_i y_{i+1} - x_{i+1} y_i)] [\iint\limits_A xy\; dx dy = \frac{1}{24} \sum_{i=0}^{n-1} (2 x_{i+1} y_{i+1} + x_{i+1} y_i + x_i y_{i+1} + 2 x_i y_i)\; (x_i y_{i+1} - x_{i+1} y_i)]

These formulas, and the terms accounting for the location of the centroid, are in the function inertia.

def inertia(pts):
  'Moments and product of inertia about centroid.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  sxx = syy = sxy = 0
  a = area(pts)
  cx, cy = centroid(pts)
  for i in range(len(pts) - 1):
    sxx += (y[i]**2 + y[i]*y[i+1] + y[i+1]**2)*(x[i]*y[i+1] - x[i+1]*y[i])
    syy += (x[i]**2 + x[i]*x[i+1] + x[i+1]**2)*(x[i]*y[i+1] - x[i+1]*y[i])
    sxy += (x[i]*y[i+1] + 2*x[i]*y[i] + 2*x[i+1]*y[i+1] + x[i+1]*y[i])*(x[i]*y[i+1] - x[i+1]*y[i])
  return sxx/12 - a*cy**2, syy/12 - a*cx**2, sxy/24 - a*cx*cy

This older post explains the use of the moment of inertia in beam bending, but I avoided the trickier bits associated with the product of inertia and principal axes. We’ll cover them in the next post.

Update Jan 23, 2018 12:42 PM  Thanks to Glenn Walker for finding an error in one of the formulas. They’re more annoying to me than mistakes in the text.

  1. Yes, this is the parallel axis theorem

  2. Yes, the product of inertia integral is definitely more complicated if you’re going to do the derivation by hand. So don’t do it by hand. Learn SymPy and you’ll be able to zip through it.