Next post Previous post

Section properties and SymPy

January 18, 2018 at 5:28 PM by Dr. Drang

I felt a little guilty about this footnote in yesterday’s post:

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.

This is entirely too much like those “it can be easily shown” tricks that math textbook writers use to avoid complicated and unintuitive manipulations. If I’m going to claim you can zip through the product of inertia, I should be able to prove it. So let’s do it.

SymPy comes with the Anaconda Python distribution, and that’s how I installed it. I believe you can get it working with Apple’s system-supplied Python, but Anaconda is so helpful in getting and maintaining a numerical/scientific Python installation, I don’t see why you’d try anything else.

If you’ve ever used a symbolic math program, like Mathematica or Maple, SymPy will seem reasonably familiar to you. My main hangup is the need in SymPy to declare certain variables as symbols before doing any other work. I understand the reason for it—SymPy needs to protect symbols from being evaluated the way regular Python variables are—but I tend to forget to declare all the symbols I need and don’t realize it until an error message appears.

That one personal quirk aside, I find SymPy easy to use for the elementary math I tend to do. The functions I use most often, like diff, integrate, expand, and factor, are easy to remember, so I don’t have to continually look things up in the documentation. And the docs are well-organized when I do have to use them.

The problem we’re going to look at is the solution of this integral for a polygonal area:

∬Axydxdy\iint\limits_A xy\; dx dy

We’ll use Green’s theorem to turn this area integral into a path integral around the polygon’s perimeter:

∬Axydxdy=∮C12x2ydy\iint\limits_A xy\; dx dy = \oint\limits_C \frac{1}{2} x^2 y\; dy

For each side of the polygon, from point (xi,yi)(x_i, y_i) to point (xi+1,yi+1)(x_{i+1}, y_{i+1}), the line segment defining the perimeter can be expressed in parametric form,

x=xi+(xi+1−xi)tx = x_i + (x_{i+1} - x_i)\;t y=yi+(yi+1−yi)ty = y_i + (y_{i+1} - y_i)\;t

which means

dy=(yi+1−yi)dtdy = (y_{i+1} - y_i)\; dt

Now we’re ready to use SymPy to evaluate and simplify the integral for a single line segment. To make the typing go faster as I used SymPy, which I ran interactively in Jupyter console session, I decided to use 0 for subscript ii and 1 for subscript i+1i+1. Here’s a transcript of the session, where I’ve broken up long lines to make it easier to read:

In [1]: from sympy import *

In [2]: x, y, x_0, x_1, y_0, y_1, t = symbols('x y x_0 x_1 y_0 y_1 t')

In [3]: x = x_0 + (x_1 - x_0)*t

In [4]: y = y_0 + (y_1 - y_0)*t

In [5]: full = integrate(x**2*y/2*diff(y, t), (t, 0, 1))

In [6]: full
Out[6]: -x_0**2*y_0**2/8 + x_0**2*y_0*y_1/12 + x_0**2*y_1**2/24
 - x_0*x_1*y_0**2/12 + x_0*x_1*y_1**2/12 - x_1**2*y_0**2/24
 - x_1**2*y_0*y_1/12 + x_1**2*y_1**2/8

In [7]: part = x_0**2*y_0*y_1/12 + x_0**2*y_1**2/24 - x_0*x_1*y_0**2/12
 + x_0*x_1*y_1**2/12 - x_1**2*y_0**2/24 - x_1**2*y_0*y_1/12

In [8]: factor(part)
Out[8]: (x_0*y_1 - x_1*y_0)*(2*x_0*y_0 + x_0*y_1 + x_1*y_0 + 2*x_1*y_1)/24

In [9]: print(latex(_))
\frac{1}{24} \left(x_{0} y_{1} - x_{1} y_{0}\right) \left(2 x_{0} y_{0}
 + x_{0} y_{1} + x_{1} y_{0} + 2 x_{1} y_{1}\right)

We start by importing everything from SymPy and defining all the symbols needed. Then we define the parametric equations of the line segment in In[3] and In[4].

In[5] does a lot of work. We define the integrand inside the integrate function and tell it to integrate that expression over tt from 0 to 1 (i.e., from (x0,y0)(x_0, y_0) to (x1,y1)(x_1, y_1)). Note that we didn’t need to explicitly enter the expressions for xx, yy, or dydy; SymPy did all the substitution for us, including the differentiation.

I called the result of the integration full because it contains every term of the integration. But we learned in the last post that the leading and trailing terms get cancelled out when we sum over all the segments of the polygon. So I copied just the inner terms from full and pasted them into In[7] to define a new expression, called part.

In[8] then factors part to get a more compact expression, and In[9] converts it to a LaTeX expression, so I can render it nicely here:

124(x0y1−x1y0)(2x0y0+x0y1+x1y0+2x1y1)\frac{1}{24} \left(x_{0} y_{1} - x_{1} y_{0}\right) \left(2 x_{0} y_{0} + x_{0} y_{1} + x_{1} y_{0} + 2 x_{1} y_{1}\right)

With a quick search-and-replace to convert the subscripts to their more general forms, we get the expression presented in the last post (albeit with the terms in a different order):

∬Axydxdy=124∑i=0n−1(xiyi+1−xi+1yi)(2xiyi+xiyi+1+xi+1yi+2xi+1yi+1)\iint\limits_A xy\; dx dy = \frac{1}{24} \sum_{i=0}^{n-1} \left(x_{i} y_{i+1} - x_{i+1} y_{i}\right) \left(2 x_{i} y_{i} + x_{i} y_{i+1} + x_{i+1} y_{i} + 2 x_{i+1} y_{i+1}\right)

SymPy didn’t do everything for us. We had to figure out the Green’s function transformation and recognize the cancellation of the leading and trailing terms of full. But it did all the boring stuff, which is its real value.

Next post Previous post