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
from scipy.integrate import nquad

Drawing Bézier curves

Reference: This YouTube video.

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 = fig.add_subplot(111, projection='3d')
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