# Calculus tools and demonstrations¶

By Max Kapur. Updated May. 27, 2020.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from time import perf_counter
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import binom


## Drawing Bézier curves

Fix a list of $n$ points in $\mathbb{R}^m$. We would like to create a curve that proceeds smoothly from point $x_1$ to $x_n$ so that when the curve is closest to $x_i$, its normal vector points toward the point $x_i$. At the endpoints $x_1$ and $x_n$, the curve should be tangent to the straight lines connecting $x_1$ and $x_2$, and $x_{n-1}$ and $x_n$, respectively. This is a Bézier curve, useful in design and illustration.

The solution is an iterative parameterization process. First, create vector functions of $t$ that make a linear interpolation between each point in the sequence (except the final point $x_n$) and the point immediately following it, yielding a vector of (vector) functions $y = \left<y_1(t), y_2(t), \dots ,y_{n-1}(t)\right>$. Apply this same procedure to the $n-1$ vectors in $y$, yielding $z = \left<z_1(t), z_2(t), \dots ,z_{n-2}(t)\right>$, and keep iterating until the resulting list contains only one vector function, $r(t)$. This is the Bézier curve, a polynomial of degree $n-1$.

In [3]:
fig, ax = plt.subplots(1,1,figsize=(6,6))
x = np.array([[1,1],[2, 1.5], [2.5, 3]])
plt.scatter(*x.T, color='black')
plt.xlim(.9, 2.7)
for i, xy in enumerate(x):
plt.annotate(r'$x_{}$'.format(i+1), xy+[0.1,0], ha='center', va='center',size=15)


For the points shown above, we have $n=3$ and $m=2$, with

\begin{align}x_1 &= \begin{bmatrix}1 \\ 1\end{bmatrix} \\ x_2 &= \begin{bmatrix}2 \\ 1.5\end{bmatrix} \\ x_3 &= \begin{bmatrix}2.5 \\ 3\end{bmatrix} \end{align}

The first iteration is as follows:

\begin{align} \text{from } x_1 \to x_2:\quad & y_1 (t) & = (1-t) x_1 + (t) x_2 \\ \text{from } x_2 \to x_3:\quad & y_2 (t) & = (1-t) x_2 + (t) x_3 \\ \end{align}

In the second and final iteration, $r(t)$ connects $y_1$ and $y_2$:

\begin{align}r(t) & = (1-t) y_1 + (t) y_2 \\ & = (1-t) \left[(1-t) x_1 + (t) x_2\right] + (t) \left[(1-t) x_2 + (t) x_3\right] \\ & = \begin{bmatrix}- 0.5 \\1.0 \end{bmatrix} t^2 + \begin{bmatrix} 2 \\ 1\end{bmatrix} t + \begin{bmatrix} 1\\1 \end{bmatrix} \\ & = \begin{bmatrix}- 0.5 t^{2} + 2 t + 1\\1.0 t^{2} + 1 t + 1\end{bmatrix} \end{align}

This is quadratic in $t$, as expected.

Here I’ve written a Python function that draws Bézier curves of a specified resolution for an arbitrary list of points.

In [4]:
def bezier(pivots, res=10):
assert len(pivots.shape)==2

# t-values at which curve will be evaluated
t = np.linspace(0,1,res)

# primitive curves (segments between points)
p=np.stack([  np.outer((1-t), pivots[i])
+ np.outer((t), pivots[i+1])
for i in range(pivots.shape[0] - 1)])

# recursively parameterize between the curves until we've done them all
while p.shape[0] >= 2:
p=np.stack([ (p[i].T * (1-t) ).T
+(p[i+1].T * (t) ).T
for i in range(p.shape[0] - 1)])
return p.reshape(p.shape[1:])

In [15]:
four = [np.array([[1,1],[2, 1.5], [2.5, 3],[1.8, 4]]),
np.array([1.3, 1.2,2.85, 1.65,1.85, 3.4,4.05, 3.6,6.05, 3.05]).reshape(5,2),
np.array([[0.9, 3.75],
[0.4, 2.85],
[0.85, 2],
[2.65, 0.75]]),
np.array([[0.55, 2.35],
[3.3, 4.15],
[3.05, 0.95]])]

colors = plt.cm.Dark2(np.linspace(0,.49,4))
fig, ax = plt.subplots(2,2,figsize=(12,12))

for i, pivots in enumerate(four):
ax[int(i/2), i%2].plot(*bezier(pivots, res=500).T, color = colors[i])
ax[int(i/2), i%2].scatter(*pivots.T, color="black")
for j, xy in enumerate(pivots):
ax[int(i/2), i%2].annotate(j, xy + [.05,-.05], ha='center', va='center')


It works in 3D, too.

In [19]:
pivots = np.array([[0.8, 1.2, 1],
[1.55, 3.8, 1],
[2.25, 1.75, 1],
[3.4, 2.65, 1],
[3.25, 2.5, 1]])

fig = plt.figure(figsize=(12,10))
ax.view_init(45, 45)

ax.scatter(*pivots.T)
ax.plot(*bezier(pivots, res=500).T, color = 'crimson')

Out[19]:
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x17fd62d3ac8>]

An explicit definition of the Bézier curve involves a binomial expansion: $$r(t) = \sum_{i=0}^n \binom{n}{i}(1-t)^{n-i}t^i x_i$$ Here is another, perhaps more legible, Python function that uses this explicit definition instead:

In [5]:
def bezier2(pivots, res=10):
assert len(pivots.shape)==2

t = np.linspace(0,1, res)
n = pivots.shape[0]-1

return sum([binom(n,i) * np.outer((1-t)**(n-i) * t**i, x)
for i, x in enumerate(pivots)])


Let’s compare the performance of the two functions.

In [69]:
np.random.seed(2020)

# problems of various sizes
test_sets = [np.random.rand(n,m) for n in range(3,20) for m in range(2,20) for i in range(10)]
len(test_sets)

Out[69]:
3060
In [70]:
sizes = []
times_1 = []
times_2 = []

for i in test_sets:
sizes.append(np.product(i.shape))

tic = perf_counter()
bezier(i, res=100)
toc = perf_counter()
times_1.append(toc-tic)

tic = perf_counter()
bezier2(i, res=100)
toc = perf_counter()
times_2.append(toc-tic)

In [73]:
sizes, times_1, times_2 = [np.array(i) for i in [sizes, times_1, times_2]]
times_1 *= 1e3
times_2 *= 1e3

z1 = np.polyfit(sizes, times_1, 1)
p1 = np.poly1d(z1)

z2 = np.polyfit(sizes, times_2, 1)
p2 = np.poly1d(z2)

plt.figure(figsize = (10,8))
plt.xlim(0,405)

plt.scatter(sizes, times_1, label = 'list comprehension', marker = '+', alpha = .7, color='tab:olive')
plt.plot(plt.xlim(), p1(plt.xlim()), color='tab:olive')

plt.scatter(sizes, times_2, label = 'binomial expansion', marker = '+', alpha = .7, color='tab:purple')
plt.plot(plt.xlim(), p2(plt.xlim()), color='tab:purple')

plt.xlabel(r'$n\times m$')
plt.ylabel('time (ms)')

plt.legend()
plt.title('Bézier curves: Two implementations', size=16)
plt.xlim()

Out[73]:
(0.0, 405.0)

Looks like the binomial expansion version is much more performant, which is to be expected given the enhancements built into scipy.special.binom() that we miss when expanding the binomial “manually.”

## Moment of inertia of a solid ellipsoid¶

References:

Consider a solid ellipsoid in $\mathbb{R}^3$ defined by

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$$

having mass $m$ and constant density $\lambda$. I derive expressions for the moments of inertia about the coordinate axes in terms of $m$, $a$, $b$, and $c$.

### Volume of ellipsoid¶

The volume of the ellipsoid is given by

$$V_E = \iiint_E \, dV$$

where $D$ is the ellipsoid itself.

Under the transformation

$$x = ua \\ y = vb \\ z = wc$$

the ellipsoid becomes a unit sphere. Then

$$V_E = \iiint_{E'} \left\vert J \right\vert \, dV$$

where $E'$ is the unit sphere. Here

$$\left\vert J \right\vert = \left\vert \begin{matrix} a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{matrix} \right\vert = abc$$

is a constant, so

\begin{align} V_E & = abc \iiint_{E'} dV \\ & = abc \cdot \text{(volume of unit sphere)} \\ & = \frac{4}{3} \pi abc \end{align}

### Moment of inertia about $z$-axis¶

The moment of inertia is

$$I_z = \int r^2 \, dm$$

where $r$ indicates the distance from the axis of rotation and $dm= \lambda dV$. Here $r^2$ (not to be confused with $r$ later used as an argument in spherical coordinates) is the squared distance from the $z$-axis, $x^2 + y^2$, so

$$I_z = \iiint_E x^2 + y^2 \, dV$$

where E is the ellipsoid. Or, applying the same transformation as above,

\begin{align} I_z &= abc \iiint_{E'} x^2 + y^2 \, dV \\ & = abc \iiint_{E'} a^2 u^2 + b^2 v^2 \, dV \end{align}

where $E'$ is the unit sphere.

Now, transform to spherical coordinates as follows:

$$\begin{bmatrix} u \\ v \\ w \end{bmatrix} = \begin{bmatrix} r \sin \theta \cos \phi \\ r \sin \theta \sin \phi \\ r \cos \theta \end{bmatrix}$$

Under the transformation, the unit sphere is becomes the cubic region $D$ as follows:

D: \begin{align} r & \in [0,1] \\ \theta & \in [0,\pi] \\ \phi & \in [0,2\pi] \end{align}

Also, for the transformation to spherical coordinates

$$\left\vert J \right\vert = \left\vert \begin{matrix} \sin \theta \cos \phi & r \cos \theta \cos \phi & -r \sin \theta\sin\phi \\ \sin \theta \sin \phi & r \cos \theta \sin \phi & -r \sin \theta\cos\phi \\ \cos \theta & -r \sin \theta & 0 \\ \end{matrix} \right\vert = r^2 \sin\theta$$

(trust me). So,

\begin{align} I_z &= abc \iiint_{D} \left( a^2 u^2 + b^2 v^2\right) \left\vert J \right\vert \, dV \\ & = abc \int_0^{2\pi} \int_0^{\pi} \int_0^1 \left( a^2 r^2 \sin^2\theta + b^2 r^2 \sin^2\theta \right) r^2 \sin \theta \, dr\,d \theta \,d\phi \\ & = a^3 bc \left(A\right) + ab^3c \left(A\right) \end{align}

by factoring out constants and defining

\begin{align} A & \equiv \int_0^{2\pi} \int_0^{\pi} \int_0^1 r^4 \sin^3\theta \cos^2 \phi \, dr\,d \theta \,d\phi &= \frac{4}{15} \pi\\ B & \equiv \int_0^{2\pi} \int_0^{\pi} \int_0^1 r^4 \sin^3\theta \sin^2 \phi \, dr\,d \theta \,d\phi &= \frac{4}{15} \pi \end{align}

These integrals can be computed explicitly, but here is a quick verification using scipy.integrate:

In [5]:
def A(r, t, p):
return r**4 * np.sin(t)**3 * np.cos(p)**2

def B(r, t, p):
return r**4 * np.sin(t)**3 * np.sin(p)**2

print(nquad(A, [[0, 1], [0, np.pi], [0, np.pi*2]])[0])
print(nquad(B, [[0, 1], [0, np.pi], [0, np.pi*2]])[0])
print(4*np.pi / 15)

0.837758040957278
0.8377580409572783
0.8377580409572781


Finally,

\begin{align} I_z &= \lambda \left( \frac{4}{15} \pi \right) \left(a^3 bc + ab^3 c \right)\\ &= \frac{\lambda}{5} \left( \frac{4}{3} \pi abc \right) \left(a^2 + b^2 \right)\\ m &=\lambda V_E = \lambda \left(\frac{4}{3} \pi abc\right) \\ & \implies\quad I_z = \frac{m}{5}\left(a^2 + b^2\right) \end{align}

### Moments of inertia about principle axes¶

By symmetry, the moments of inertia about the other axes can be obtained by exchanging dimensions:

\begin{align} I_x & = \frac{m}{5}\left(b^2 + c^2\right) \\ I_y & = \frac{m}{5}\left(a^2 + c^2\right) \\ I_z & = \frac{m}{5}\left(a^2 + b^2\right) \end{align}
In [9]:
!jupyter nbconvert calc.ipynb

[NbConvertApp] Converting notebook calc.ipynb to html
[NbConvertApp] Writing 454040 bytes to calc.html