The mass and centre of mass of a tetrahedron

This example finds the mass and centre of mass of the tetrahedron bounded by the coordinate axes and the plane $x+y+z=1$ with density $\rho = \rho(x, y, z)$ where $\rho(x,y,z)$ is provided as a lambda function. We test it with the functions $\rho = 1$, $\rho = x$ and $\rho = x^2+y^2+z^2$.

The mass may be written as a triple integral of the density over the volume of the tetrahedron: $$ m = \int_V \rho(x, y, z)\;\mathrm{d}V = \int_0^1\int_0^{1-x}\int_0^{1-x-y} \rho(x, y, z) \;\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x, $$ and the coordinates of the centre of mass are given by $$ m\bar{x} = \int_V x\rho(x,y,z)\;\mathrm{d}V, \quad m\bar{y} = \int_V y\rho(x,y,z)\;\mathrm{d}V, \quad m\bar{z} = \int_V z\rho(x,y,z)\;\mathrm{d}V. $$

The following program uses scipy.integrate.tplquad to perform the necessary integrations (which can also be solved analytically).

import numpy as np
from scipy.integrate import tplquad

# The integration limits on x, y, z:
a, b = 0, 1
gfun, hfun = lambda x: 0, lambda x: 1 - x
qfun, rfun = lambda x, y: 0, lambda x, y: 1 - x - y
lims = (a, b, gfun, hfun, qfun, rfun)

# The three different density functions
rhos = [lambda x, y, z: 1,
        lambda x, y, z: x,
        lambda x, y, z: x**2 + y**2 + z**2]

for rho in rhos:
    # The mass as a triple integral of rho over the volume
    m, _ = tplquad(rho, *lims)
    # The centre of mass (xbar, ybar, zbar)
    mxbar, _ = tplquad(lambda x, y, z: x * rho(x,y,z), *lims)
    mybar, _ = tplquad(lambda x, y, z: y * rho(x,y,z), *lims)
    mzbar, _ = tplquad(lambda x, y, z: z * rho(x,y,z), *lims)
    xbar, ybar, zbar = mxbar / m, mybar / m, mzbar / m

    print('mass = {:g}, CofM = ({:g}, {:g}, {:g})'.format(m, xbar, ybar, zbar))

Note that the six arguments representing the limits on the triple integral (two constants and two pairs of lambda functions) have been packed into a tuple, lims (the parentheses are optional here).

The output is:

mass = 0.166667, CofM = (0.25, 0.25, 0.25)
mass = 0.0416667, CofM = (0.4, 0.2, 0.2)
mass = 0.05, CofM = (0.277778, 0.277778, 0.277778)