chapter7b.tex

 \chapter{Applications of Tensor}


\section{Common Definitions in Tensor Notation\label{subCommonDefinitions}}

 The trace of a matrix $\vector{A}$ representing a rank-2
tensor is:
\begin{equation}
\mathrm{tr}\left(\vector{A}\right)=A_{ii}
\end{equation}


 For a $3\times3$ matrix representing a rank-2 tensor in
a 3D space, the determinant is:
\begin{equation}
\mathrm{det}\left(\vector{A}\right)=\begin{vmatrix}\begin{array}{ccc}
A_{11} & A_{12} & A_{13}\\
A_{21} & A_{22} & A_{23}\\
A_{31} & A_{32} & A_{33}
\end{array}\end{vmatrix}=\epsilon_{ijk}A_{1i}A_{2j}A_{3k}=\epsilon_{ijk}A_{i1}A_
{j2}A_{k3}
\end{equation}
where the last two equalities represent the expansion of the determinant
by row and by column. Alternatively
\begin{equation}
\mathrm{det}\left(\vector{A}\right)=\dfrac{1}{3!}\epsilon_{ijk}\epsilon_{lmn}A_{
il}A_{jm}A_{kn}
\end{equation}


 For an $n\times n$ matrix representing a rank-2 tensor
in $\bbR^n$, the determinant is:
\begin{equation}
\mathrm{det}\left(\vector{A}\right)=\epsilon_{i_{1}\cdots i_{n}}A_{1i_{1}}\ldots 
A_{ni_{n}}=\epsilon_{i_{1}\cdots i_{n}}A_{i_{1}1}\ldots 
A_{i_{n}n}=\dfrac{1}{n!}\epsilon_{i_{1}\cdots i_{n}}\,\epsilon_{j_{1}\cdots 
j_{n}}A_{i_{1}j_{1}}\ldots A_{i_{n}j_{n}}
\end{equation}


 The inverse of a matrix $\vector{A}$ representing a rank-2
tensor is:
\begin{equation}
\left[\vector{A}^{-1}\right]_{ij}=\dfrac{1}{2\,\mathrm{det}\left(\vector{A}
\right)}\epsilon_{jmn}\,\epsilon_{ipq}A_{mp}A_{nq}
\end{equation}


 The multiplication of a matrix $\vector{A}$ by a vector
$\vector{b}$ as defined in linear algebra is:
\begin{equation}
\left[\vector{A}\vector{b}\right]_{i}=A_{ij}b_{j}
\end{equation}
It should be noticed that here we are using matrix notation. The multiplication
operation, according to the symbolic notation of tensors, should be
denoted by a dot between the tensor and the vector, i.e. 
$\vector{A}\vector{\cdot b}$.\footnote{The matrix multiplication in matrix 
notation is equivalent to a dot
product operation in tensor notation.}

 The multiplication of two $n\times n$ matrices $\vector{A}$
and $\vector{B}$ as defined in linear algebra is:
\begin{equation}
\left[\vector{A}\vector{B}\right]_{ik}=A_{ij}B_{jk}\label{eqMatrixMultiplication
}
\end{equation}
Again, here we are using matrix notation; otherwise a dot should be
inserted between the two matrices.

 The dot product of two vectors is:
\begin{equation}
\vector{A}\cdot\vector{B=}\delta_{ij}A_{i}B_{j}=A_{i}B_{i}\label{eqDotProduct}
\end{equation}
The readers are referred to $\S$ \ref{secInnerProduct} for a more
general definition of this type of product that includes higher rank
tensors.

 The cross product of two vectors is:
\begin{equation}
\left[\vector{A}\times\vector{B}\right]_{i}=\epsilon_{ijk}A_{j}B_{k}\label{
EqCrossProduct}
\end{equation}


 The scalar triple product of three vectors is:
\begin{equation}
\vector{A}\cdot\left(\vector{B}\times\vector{C}\right)=\begin{vmatrix}
\begin{array}{ccc}
A_{1} & A_{2} & A_{3}\\
B_{1} & B_{2} & B_{3}\\
C_{1} & C_{2} & C_{3}
\end{array}\end{vmatrix}=\epsilon_{ijk}A_{i}B_{j}C_{k}\label{
EqScalarTripleProduct}
\end{equation}


 The vector triple product of three vectors is:
\begin{equation}
\left[\vector{A}\times\left(\vector{B}\times\vector{C}\right)\right]_{i}
=\epsilon_{ijk}\epsilon_{klm}A_{j}B_{l}C_{m}
\end{equation}
% 
% 
% 
% \section{Scalar Invariants of Tensors}
% 
%  In the following we list and write in tensor notation a
% number of invariants of low rank tensors which have special importance
% due to their widespread applications in vector and tensor calculus.
% All These invariants are scalars.
% 
%  The value of a scalar (rank-0 tensor), which consists of
% a magnitude and a sign, is invariant under coordinate transformation.
% 
%  An invariant of a vector (rank-1 tensor) under coordinate
% transformations is its magnitude, i.e. length (the direction is also
% invariant but it is not scalar!).\footnote{In fact the magnitude alone is 
% invariant under coordinate transformations
% even for pseudo vectors because it is a scalar.}
% 
%  The main three independent scalar invariants of a rank-2
% tensor $\vector{A}$ under change of basis are:
% \begin{equation}
% I=\mathrm{tr}\left(\vector{A}\right)=A_{ii}
% \end{equation}
% \begin{equation}
% II=\mathrm{tr}\left(\vector{A}^{2}\right)=A_{ij}A_{ji}
% \end{equation}
% \begin{equation}
% III=\mathrm{tr}\left(\vector{A}^{3}\right)=A_{ij}A_{jk}A_{ki}
% \end{equation}
% 
% 
%  Different forms of the three invariants of a rank-2 tensor
% $\vector{A}$, which are also widely used, are:
% \begin{equation}
% I_{1}=I=A_{ii}
% \end{equation}
% \begin{equation}
% I_{2}=\dfrac{1}{2}\left(I^{2}-II\right)=\dfrac{1}{2}\left(A_{ii}A_{jj}-A_{ij}A_{ji
% }\right)
% \end{equation}
% \begin{equation}
% I_{3}=\mathrm{det}\left(\vector{A}\right)=\dfrac{1}{3!}\left(I^{3}-3I\,\,
% II+2III\right)=\dfrac{1}{3!}\epsilon_{ijk}\epsilon_{pqr}A_{ip}A_{jq}A_{kr}
% \end{equation}
% 
% 
%  The invariants $I$, $II$ and $III$ can similarly be
% defined in terms of the invariants $I_{1}$, $I_{2}$ and $I_{3}$
% as follow:
% \begin{equation}
% I=I_{1}
% \end{equation}
% \begin{equation}
% II=I_{1}^{2}-2I_{2}
% \end{equation}
% \begin{equation}
% III=I_{1}^{3}-3I_{1}I_{2}+3I_{3}
% \end{equation}
% 
% 
%  Since the determinant of a matrix representing a rank-2
% tensor is invariant, then if the determinant vanishes in one coordinate
% system it will vanish in all coordinate systems and vice versa. Consequently,
% if a rank-2 tensor is invertible in a particular coordinate system,
% it is invertible in all coordinate systems.
% 
%  Ten joint invariants between two rank-2 tensors, $\vector{A}$
% and $\vector{B}$, can be formed; these are: 
% $\mathrm{tr}\left(\vector{A}\right)$,
% $\mathrm{tr}\left(\vector{B}\right)$, $\mathrm{tr}\left(\vector{A}^{2}\right)$,
% $\mathrm{tr}\left(\vector{B}^{2}\right)$, 
% $\mathrm{tr}\left(\vector{A}^{3}\right)$,
% $\mathrm{tr}\left(\vector{B}^{3}\right)$, 
% $\mathrm{tr}\left(\vector{A}\cdot\vector{B}\right)$,
% $\mathrm{tr}\left(\vector{A}^{2}\cdot\vector{B}\right)$, 
% $\mathrm{tr}\left(\vector{A}\cdot\vector{B}^{2}\right)$
% and $\mathrm{tr}\left(\vector{A}^{2}\cdot\vector{B}^{2}\right)$.


\section{Common Differential Operations in Tensor Notation}

 Here we present the most common differential operations
as defined by tensor notation. These operations are mostly based on
the various types of interaction between the vector differential operator
nabla $\nabla$ with tensors of different ranks as well as interaction
with other types of operation like dot and cross products.

The operator $\nabla$ is essentially a spatial partial differential
operator defined in Cartesian coordinate systems by:
\begin{equation}
\nabla_{i}=\dfrac{\partial}{\partial x_{i}}
\end{equation}


 The gradient of a differentiable scalar function of position
$f$ is a vector given by:
\begin{equation}
\left[\nabla f\right]_{i}=\nabla_{i}f=\dfrac{\partial f}{\partial 
x_{i}}=\partial_{i}f=f_{,i}\label{eqGrad}
\end{equation}


 The gradient of a differentiable vector function of position
$\vector{A}$ (which is the outer product, as defined in $\S$ 
\ref{subTensorMultiplication},
between the $\nabla$ operator and the vector) is a rank-2 tensor
defined by:
\begin{equation}
\left[\nabla\vector{A}\right]_{ij}=\partial_{i}A_{j}\label{eqGrad2}
\end{equation}


 The gradient operation is distributive but not commutative
or associative:
\begin{equation}
\nabla\left(f+h\right)=\nabla f+\nabla h
\end{equation}
\begin{equation}
\nabla f\ne f\nabla
\end{equation}
\begin{equation}
\left(\nabla f\right)h\ne\nabla\left(fh\right)
\end{equation}
where $f$ and $h$ are differentiable scalar functions of position.

 The divergence of a differentiable vector $\vector{A}$
is a scalar given by:
\begin{equation}
\nabla\cdot\vector{A}=\delta_{ij}\dfrac{\partial A_{i}}{\partial 
x_{j}}=\dfrac{\partial A_{i}}{\partial 
x_{i}}=\nabla_{i}A_{i}=\partial_{i}A_{i}=A_{i,i}\label{eqDiv}
\end{equation}
The divergence operation can also be viewed as taking the gradient
of the vector followed by a contraction. Hence, the divergence of
a vector is invariant because it is the trace of a rank-2 tensor.\footnote{It 
may also be argued that the divergence of a vector is a scalar
and hence it is invariant.}

 The divergence of a differentiable rank-2 tensor $\vector{A}$
is a vector defined in one of its forms by:
\begin{equation}
\left[\nabla\cdot\vector{A}\right]_{i}=\partial_{j}A_{ji}\label{eqDiv2}
\end{equation}
and in another form by
\begin{equation}
\left[\nabla\cdot\vector{A}\right]_{j}=\partial_{i}A_{ji}\label{eqDiv3}
\end{equation}
These two different forms can be given, respectively, in symbolic
notation by:
\begin{equation}
\nabla\cdot\vector{A}\,\,\,\,\,\,\,\,\,\,\&\,\,\,\,\,\,\,\,\,\,
\nabla\cdot\vector{A}^{T}
\end{equation}
where $\vector{A}^{T}$ is the transpose of $\vector{A}$. More generally,
the divergence of a tensor of rank $n\ge2$, which is a tensor of
rank-($n-1$), can be defined in several forms, which are different
in general, depending on the combination of the contracted indices.

 The divergence operation is distributive but not commutative
or associative:
\begin{equation}
\nabla\cdot\left(\vector{A}+\vector{B}\right)=\nabla\cdot\vector{A}
+\nabla\cdot\vector{B}
\end{equation}
\begin{equation}
\nabla\cdot\vector{A}\ne\vector{A}\cdot\nabla
\end{equation}
\begin{equation}
\nabla\cdot\left(f\vector{A}\right)\ne\nabla f\cdot\vector{A}
\end{equation}
where $\vector{A}$ and $\vector{B}$ are differentiable tensor functions
of position.

 The curl of a differentiable vector $\vector{A}$ is a
vector given by:
\begin{equation}
\left[\nabla\times\vector{A}\right]_{i}=\epsilon_{ijk}\dfrac{\partial 
A_{k}}{\partial 
x_{j}}=\epsilon_{ijk}\nabla_{j}A_{k}=\epsilon_{ijk}\partial_{j}A_{k}=\epsilon_{
ijk}A_{k,j}\label{EqCurl}
\end{equation}


 The curl operation may be generalized to tensors of rank
$>1$, and hence the curl of a differentiable rank-2 tensor $\vector{A}$
can be defined as a rank-2 tensor given by:
\begin{equation}
\left[\nabla\times\vector{A}\right]_{ij}=\epsilon_{imn}\partial_{m}A_{nj}
\end{equation}


 The curl operation is distributive but not commutative
or associative:
\begin{equation}
\nabla\times\left(\vector{A}+\vector{B}\right)=\nabla\times\vector{A}
+\nabla\times\vector{B}
\end{equation}
\begin{equation}
\nabla\times\vector{A}\ne\vector{A}\times\nabla
\end{equation}
\begin{equation}
\nabla\times\left(\vector{A\times 
B}\right)\ne\left(\nabla\times\vector{A}\right)\times\vector{B}
\end{equation}


 The Laplacian scalar operator, also called the harmonic
operator, acting on a differentiable scalar $f$ is given by:
\begin{equation}
\Delta f=\nabla^{2}f=\delta_{ij}\dfrac{\partial^{2}f}{\partial x_{i}\partial 
x_{j}}=\dfrac{\partial^{2}f}{\partial x_{i}\partial 
x_{i}}=\nabla_{ii}f=\partial_{ii}f=f_{,ii}\label{eqLaplacian}
\end{equation}


 The Laplacian operator acting on a differentiable vector
$\vector{A}$ is defined for each component of the vector similar
to the definition of the Laplacian acting on a scalar, that is
\begin{equation}
\left[\nabla^{2}\vector{A}\right]_{i}=\partial_{jj}A_{i}\label{eqLaplacian2}
\end{equation}


 The following scalar differential operator is commonly
used in science (e.g. in fluid dynamics):
\begin{equation}
\vector{A}\cdot\nabla=A_{i}\nabla_{i}=A_{i}\dfrac{\partial}{\partial 
x_{i}}=A_{i}\partial_{i}\label{eqANabla}
\end{equation}
where $\vector{A}$ is a vector. As indicated earlier, the order of
$A_{i}$ and $\partial_{i}$ should be respected.

 The following vector differential operator also has common
applications in science:
\begin{equation}
\left[\vector{A}\times\nabla\right]_{i}=\epsilon_{ijk}A_{j}\partial_{k}
\end{equation}


 The differentiation of a tensor increases its rank by one,
by introducing an extra covariant index, unless it implies a contraction
in which case it reduces the rank by one. Therefore the gradient of
a scalar is a vector and the gradient of a vector is a rank-2 tensor
($\partial_{i}A_{j}$), while the divergence of a vector is a scalar
and the divergence of a rank-2 tensor is a vector ($\partial_{j}A_{ji}$
or $\partial_{i}A_{ji}$). This may be justified by the fact that
$\nabla$ is a vector operator. On the other hand the Laplacian operator
does not change the rank since it is a scalar operator; hence the
Laplacian of a scalar is a scalar and the Laplacian of a vector is
a vector.





\section{Common Identities in Vector and Tensor Notation}

 Here we present some of the widely used identities of vector
calculus in the traditional vector notation and in its equivalent
tensor notation. In the following bullet points, $f$ and $h$ are
differentiable scalar fields; $\vector{A}$, $\vector{B}$, $\vector{C}$
and $\vector{D}$ are differentiable vector fields; and 
$\vector{r}=x_{i}\vector{e}_{i}$
is the position vector.


\begin{eqnarray}
\nabla\cdot\vector{r} & = & n\nonumber \\
 & \Updownarrow\\
\partial_{i}x_{i} & = & n\nonumber
\end{eqnarray}
where $n$ is the space dimension.


\begin{eqnarray}
\nabla\times\vector{r} & = & \vector{0}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\partial_{j}x_{k} & = & 0\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\left(\vector{a}\cdot\vector{r}\right) & = & \vector{a}\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(a_{j}x_{j}\right) & = & a_{i}\nonumber
\end{eqnarray}
where $\vector{a}$ is a constant vector.


\begin{eqnarray}
\nabla\cdot\left(\nabla f\right) & = & \nabla^{2}f\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(\partial_{i}f\right) & = & \partial_{ii}f\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\cdot\left(\nabla\times\vector{A}\right) & = & 0\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\partial_{i}\partial_{j}A_{k} & = & 0\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\times\left(\nabla f\right) & = & \vector{0}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\partial_{j}\partial_{k}f & = & 0\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\left(fh\right) & = & f\nabla h+h\nabla f\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(fh\right) & = & f\partial_{i}h+h\partial_{i}f\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\cdot\left(f\vector{A}\right) & = & 
f\nabla\cdot\vector{A}+\vector{A}\cdot\nabla f\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(fA_{i}\right) & = & 
f\partial_{i}A_{i}+A_{i}\partial_{i}f\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\times\left(f\vector{A}\right) & = & f\nabla\times\vector{A}+\nabla 
f\times\vector{A}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\partial_{j}\left(fA_{k}\right) & = & 
f\epsilon_{ijk}\partial_{j}A_{k}+\epsilon_{ijk}\left(\partial_{j}f\right)A_{k}
\nonumber
\end{eqnarray}



\begin{alignat}{3}
\vector{A}\cdot\left(\vector{B}\times\vector{C}\right) & = & 
\,\vector{C}\cdot\left(\vector{A}\times\vector{B}\right) & = & 
\,\vector{B}\cdot\left(\vector{C}\times\vector{A}\right)\nonumber \\
 & \Updownarrow &  & \Updownarrow\\
\epsilon_{ijk}A_{i}B_{j}C_{k} & = & \epsilon_{kij}C_{k}A_{i}B_{j} & = & 
\epsilon_{jki}B_{j}C_{k}A_{i}\nonumber
\end{alignat}



\begin{eqnarray}
\vector{A}\times\left(\vector{B}\times\vector{C}\right) & = & 
\vector{B}\left(\vector{A}\cdot\vector{C}\right)-\vector{C}\left(\vector{A}
\cdot\vector{B}\right)\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}A_{j}\epsilon_{klm}B_{l}C_{m} & = & 
B_{i}\left(A_{m}C_{m}\right)-C_{i}\left(A_{l}B_{l}\right)\nonumber
\end{eqnarray}



\begin{eqnarray}
\vector{A}\times\left(\nabla\times\vector{B}\right) & = & 
\left(\nabla\vector{B}\right)\cdot\vector{A}-\vector{A}\cdot\nabla\vector{B}
\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\epsilon_{klm}A_{j}\partial_{l}B_{m} & = & 
\left(\partial_{i}B_{m}\right)A_{m}-A_{l}\left(\partial_{l}B_{i}\right)\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\times\left(\nabla\times\vector{A}\right) & = & 
\nabla\left(\nabla\cdot\vector{A}\right)-\nabla^{2}\vector{A}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\epsilon_{klm}\partial_{j}\partial_{l}A_{m} & = & 
\partial_{i}\left(\partial_{m}A_{m}\right)-\partial_{ll}A_{i}\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\left(\vector{A}\cdot\vector{B}\right) & = & 
\vector{A}\times\left(\nabla\times\vector{B}\right)+\vector{B}
\times\left(\nabla\times\vector{A}\right)+\left(\vector{A}
\cdot\nabla\right)\vector{B}+\left(\vector{B}\cdot\nabla\right)\vector{A}
\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(A_{m}B_{m}\right) & = & 
\epsilon_{ijk}A_{j}\left(\epsilon_{klm}\partial_{l}B_{m}\right)+\epsilon_{ijk}B_
{j}\left(\epsilon_{klm}\partial_{l}A_{m}\right)+\left(A_{l}\partial_{l}\right)B_
{i}+\left(B_{l}\partial_{l}\right)A_{i}\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\cdot\left(\vector{A}\times\vector{B}\right) & = & 
\vector{B}\cdot\left(\nabla\times\vector{A}\right)-\vector{A}
\cdot\left(\nabla\times\vector{B}\right)\nonumber \\
 & \Updownarrow\\
\partial_{i}\left(\epsilon_{ijk}A_{j}B_{k}\right) & = & 
B_{k}\left(\epsilon_{kij}\partial_{i}A_{j}\right)-A_{j}\left(\epsilon_{jik}
\partial_{i}B_{k}\right)\nonumber
\end{eqnarray}



\begin{eqnarray}
\nabla\times\left(\vector{A}\times\vector{B}\right) & = & 
\left(\vector{B}\cdot\nabla\right)\vector{A}+\left(\nabla\cdot\vector{B}
\right)\vector{A}-\left(\nabla\cdot\vector{A}\right)\vector{B}-\left(\vector{A}
\cdot\nabla\right)\vector{B}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\epsilon_{klm}\partial_{j}\left(A_{l}B_{m}\right) & = & 
\left(B_{m}\partial_{m}\right)A_{i}+\left(\partial_{m}B_{m}\right)A_{i}
-\left(\partial_{j}A_{j}\right)B_{i}-\left(A_{j}\partial_{j}\right)B_{i}
\nonumber
\end{eqnarray}



\begin{eqnarray}
\left(\vector{A}\times\vector{B}\right)\cdot\left(\vector{C}\times\vector{D}
\right) & = & \begin{vmatrix}\begin{array}{cc}
\vector{A}\cdot\vector{C} & \vector{A}\cdot\vector{D}\\
\vector{B}\cdot\vector{C} & \vector{B}\cdot\vector{D}
\end{array}\end{vmatrix}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}A_{j}B_{k}\epsilon_{ilm}C_{l}D_{m} & = & 
\left(A_{l}C_{l}\right)\left(B_{m}D_{m}\right)-\left(A_{m}D_{m}\right)\left(B_{l
}C_{l}\right)\nonumber
\end{eqnarray}



\begin{eqnarray}
\left(\vector{A}\times\vector{B}\right)\times\left(\vector{C}\times\vector{D}
\right) & = & 
\left[\vector{D}\cdot\left(\vector{A}\times\vector{B}\right)\right]\vector{C}
-\left[\vector{C}\cdot\left(\vector{A}\times\vector{B}\right)\right]\vector{D}
\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\epsilon_{jmn}A_{m}B_{n}\epsilon_{kpq}C_{p}D_{q} & = & 
\left(\epsilon_{qmn}D_{q}A_{m}B_{n}\right)C_{i}-\left(\epsilon_{pmn}C_{p}A_{m}B_
{n}\right)D_{i}\nonumber
\end{eqnarray}


 In vector and tensor notations, the condition for a vector
field $\vector{A}$ to be solenoidal is:
\begin{eqnarray}
\nabla\cdot\vector{A} & = & 0\nonumber \\
 & \Updownarrow\\
\partial_{i}A_{i} & = & 0\nonumber
\end{eqnarray}


 In vector and tensor notations, the condition for a vector
field $\vector{A}$ to be irrotational is:
\begin{eqnarray}
\nabla\times\vector{A} & = & \vector{0}\nonumber \\
 & \Updownarrow\\
\epsilon_{ijk}\partial_{j}A_{k} & = & 0\nonumber
\end{eqnarray}



\section{Integral Theorems in Tensor Notation}

 The divergence theorem for a differentiable vector field
$\vector{A}$ in vector and tensor notation is:
\begin{eqnarray}
\iiint_{V}\nabla\cdot\vector{A}\,d\tau & = & 
\iint_{S}\vector{A}\cdot\vector{n}\,d\sigma\nonumber \\
 & \Updownarrow\\
\dint_{V}\partial_{i}A_{i}d\tau & = & \dint_{S}A_{i}n_{i}d\sigma\nonumber
\end{eqnarray}
where $V$ is a bounded region in an $n$D space enclosed by a generalized
surface $S$, $d\tau$ and $d\sigma$ are generalized volume and surface
elements respectively, $\vector{n}$ and $n_{i}$ are unit normal
to the surface and its $i^{th}$ component respectively, and the index
$i$ ranges over $1,\ldots,n$.

 The divergence theorem for a differentiable rank-2 tensor
field $\vector{A}$ in tensor notation for the first index is given
by:
\begin{equation}
\dint_{V}\partial_{i}A_{il}d\tau=\dint_{S}A_{il}n_{i}d\sigma
\end{equation}


 The divergence theorem for differentiable tensor fields
of higher ranks $\vector{A}$ in tensor notation for the index $k$
is:
\begin{equation}
\dint_{V}\partial_{k}A_{ij\ldots k\ldots m}d\tau=\dint_{S}A_{ij\ldots k\ldots 
m}n_{k}d\sigma
\end{equation}


 Stokes theorem for a differentiable vector field $\vector{A}$
in vector and tensor notation is:
\begin{eqnarray}
\iint_{S}\left(\nabla\times\vector{A}\right)\cdot\vector{n}\,d\sigma & = & 
\dint_{C}\vector{A}\cdot d\vector{r}\nonumber \\
 & \Updownarrow\\
\dint_{S}\epsilon_{ijk}\partial_{j}A_{k}n_{i}d\sigma & = & 
\dint_{C}A_{i}dx_{i}\nonumber
\end{eqnarray}
where $C$ stands for the perimeter of the surface $S$ and $d\vector{r}$
is the vector element tangent to the perimeter.

 Stokes theorem for a differentiable rank-2 tensor field
$\vector{A}$ in tensor notation for the first index is:
\begin{equation}
\dint_{S}\epsilon_{ijk}\partial_{j}A_{kl}n_{i}d\sigma=\dint_{C}A_{il}dx_{i}
\end{equation}


 Stokes theorem for differentiable tensor fields of higher
ranks $\vector{A}$ in tensor notation for the index $k$ is:
\begin{equation}
\dint_{S}\epsilon_{ijk}\partial_{j}A_{lm\ldots k\ldots 
n}n_{i}d\sigma=\dint_{C}A_{lm\ldots k\ldots n}dx_{k}
\end{equation}



\section{Examples of Using Tensor Techniques to Prove 
Identities\label{subProvingIdentities}}

 $\nabla\cdot\vector{r}=n$:
\begin{equation}
\begin{aligned}\nabla\cdot\vector{r} & =\partial_{i}x_{i} & 
\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{eqDiv})}\\
 & =\delta_{ii} &  & \text{(Eq. \ref{eqdxn})}\\
 & =n &  & \text{(Eq. \ref{eqdxn})}
\end{aligned}
\end{equation}


 $\nabla\times\vector{r}=\vector{0}$:
\begin{equation}
\begin{aligned}\left[\nabla\times\vector{r}\right]_{i} & 
=\epsilon_{ijk}\partial_{j}x_{k} & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & 
\text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\delta_{kj} &  & \text{(Eq. \ref{eqdxdelta})}\\
 & =\epsilon_{ijj} &  & \text{(Eq. \ref{EqIndexReplace})}\\
 & =0 &  & \text{(Eq. \ref{eqEpsilon3Definition})}
\end{aligned}
\end{equation}
Since $i$ is a free index the identity is proved for all components.

 $\nabla\left(\vector{a}\cdot\vector{r}\right)=\vector{a}$:
\begin{equation}
\begin{aligned}\left[\nabla\left(\vector{a}\cdot\vector{r}\right)\right]_{i} & 
=\partial_{i}\left(a_{j}x_{j}\right) & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, 
& \text{(Eqs. \ref{eqGrad} \& \ref{eqDotProduct})}\\
 & =a_{j}\partial_{i}x_{j}+x_{j}\partial_{i}a_{j} &  & \text{(product rule)}\\
 & =a_{j}\partial_{i}x_{j} &  & \text{(\ensuremath{a_{j}} is constant)}\\
 & =a_{j}\delta_{ji} &  & \text{(Eq. \ref{eqdxdelta})}\\
 & =a_{i} &  & \text{(Eq. \ref{EqIndexReplace})}\\
 & =\left[\vector{a}\right]_{i} &  & \text{(definition of index)}
\end{aligned}
\end{equation}
Since $i$ is a free index the identity is proved for all components.

 $\nabla\cdot\left(\nabla f\right)=\nabla^{2}f$:
\begin{equation}
\begin{aligned}\nabla\cdot\left(\nabla f\right) & =\partial_{i}\left[\nabla 
f\right]_{i} & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{eqDiv})}\\
 & =\partial_{i}\left(\partial_{i}f\right) &  & \text{(Eq. \ref{eqGrad})}\\
 & =\partial_{i}\partial_{i}f &  & \text{(rules of differentiation)}\\
 & =\partial_{ii}f &  & \text{(definition of 2nd derivative)}\\
 & =\nabla^{2}f &  & \text{(Eq. \ref{eqLaplacian})}
\end{aligned}
\end{equation}


 $\nabla\cdot\left(\nabla\times\vector{A}\right)=0$:
\begin{equation}
\begin{aligned}\nabla\cdot\left(\nabla\times\vector{A}\right) & 
=\partial_{i}\left[\nabla\times\vector{A}\right]_{i} & 
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{eqDiv})}\\
 & =\partial_{i}\left(\epsilon_{ijk}\partial_{j}A_{k}\right) &  & \text{(Eq. 
\ref{EqCurl})}\\
 & =\epsilon_{ijk}\partial_{i}\partial_{j}A_{k} &  & 
\text{(\ensuremath{\partial}\ not acting on \ensuremath{\epsilon})}\\
 & =\epsilon_{ijk}\partial_{j}\partial_{i}A_{k} &  & \text{(continuity 
condition)}\\
 & =-\epsilon_{jik}\partial_{j}\partial_{i}A_{k} &  & \text{(Eq. 
\ref{EqEpsilonCycle})}\\
 & =-\epsilon_{ijk}\partial_{i}\partial_{j}A_{k} &  & \text{(relabeling dummy 
indices \ensuremath{i} and \ensuremath{j})}\\
 & =0 &  & \text{(since 
\ensuremath{\epsilon_{ijk}\partial_{i}\partial_{j}A_{k}=-\epsilon_{ijk}\partial_
{i}\partial_{j}A_{k}})}
\end{aligned}
\end{equation}
This can also be concluded from line three by arguing that: since
by the continuity condition $\partial_{i}$ and $\partial_{j}$ can
change their order with no change in the value of the term while a
corresponding change of the order of $i$ and $j$ in $\epsilon_{ijk}$
results in a sign change, we see that each term in the sum has its
own negative and hence the terms add up to zero (see Eq. 
\ref{eqPermutingTwoFactors}).

 $\nabla\times\left(\nabla f\right)=\vector{0}$:
\begin{equation}
\begin{aligned}\left[\nabla\times\left(\nabla f\right)\right]_{i} & 
=\epsilon_{ijk}\partial_{j}\left[\nabla f\right]_{k} & 
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\partial_{j}\left(\partial_{k}f\right) &  & \text{(Eq. 
\ref{eqGrad})}\\
 & =\epsilon_{ijk}\partial_{j}\partial_{k}f &  & \text{(rules of 
differentiation)}\\
 & =\epsilon_{ijk}\partial_{k}\partial_{j}f &  & \text{(continuity condition)}\\
 & =-\epsilon_{ikj}\partial_{k}\partial_{j}f &  & \text{(Eq. 
\ref{EqEpsilonCycle})}\\
 & =-\epsilon_{ijk}\partial_{j}\partial_{k}f &  & \text{(relabeling dummy 
indices \ensuremath{j} and \ensuremath{k})}\\
 & =0 &  & \text{(since 
\ensuremath{\epsilon_{ijk}\partial_{j}\partial_{k}f=-\epsilon_{ijk}\partial_{j}
\partial_{k}f})}
\end{aligned}
\end{equation}
This can also be concluded from line three by a similar argument to
the one given in the previous point. Because $\left[\nabla\times\left(\nabla 
f\right)\right]_{i}$
is an arbitrary component, then each component is zero.

 $\nabla\left(fh\right)=f\nabla h+h\nabla f$:
\begin{equation}
\begin{aligned}\left[\nabla\left(fh\right)\right]_{i} & 
=\partial_{i}\left(fh\right) & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. 
\ref{eqGrad})}\\
 & =f\partial_{i}h+h\partial_{i}f &  & \text{(product rule)}\\
 & =\left[f\nabla h\right]_{i}+\left[h\nabla f\right]_{i} &  & \text{(Eq. 
\ref{eqGrad})}\\
 & =\left[f\nabla h+h\nabla f\right]_{i} &  & \text{(Eq. 
\ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.

 
$\nabla\cdot\left(f\vector{A}\right)=f\nabla\cdot\vector{A}+\vector{A}
\cdot\nabla f$:
\begin{equation}
\begin{aligned}\nabla\cdot\left(f\vector{A}\right) & 
=\partial_{i}\left[f\vector{A}\right]_{i} & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & 
\text{(Eq. \ref{eqDiv})}\\
 & =\partial_{i}\left(fA_{i}\right) &  & \text{(definition of index)}\\
 & =f\partial_{i}A_{i}+A_{i}\partial_{i}f &  & \text{(product rule)}\\
 & =f\nabla\cdot\vector{A}+\vector{A}\cdot\nabla f &  & \text{(Eqs. \ref{eqDiv} 
\& \ref{eqANabla})}
\end{aligned}
\end{equation}


 $\nabla\times\left(f\vector{A}\right)=f\nabla\times\vector{A}+\nabla 
f\times\vector{A}$:
\begin{equation}
\begin{aligned}\left[\nabla\times\left(f\vector{A}\right)\right]_{i} & 
=\epsilon_{ijk}\partial_{j}\left[f\vector{A}\right]_{k} & 
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\partial_{j}\left(fA_{k}\right) &  & \text{(definition of 
index)}\\
 & 
=f\epsilon_{ijk}\partial_{j}A_{k}+\epsilon_{ijk}\left(\partial_{j}f\right)A_{k} 
&  & \text{(product rule \& commutativity)}\\
 & =f\epsilon_{ijk}\partial_{j}A_{k}+\epsilon_{ijk}\left[\nabla 
f\right]_{j}A_{k} &  & \text{(Eq. \ref{eqGrad})}\\
 & =\left[f\nabla\times\vector{A}\right]_{i}+\left[\nabla 
f\times\vector{A}\right]_{i} &  & \text{(Eqs. \ref{EqCurl} \& 
\ref{EqCrossProduct})}\\
 & =\left[f\nabla\times\vector{A}+\nabla f\times\vector{A}\right]_{i} &  & 
\text{(Eq. \ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.

 
$\vector{A}\cdot\left(\vector{B}\times\vector{C}\right)=\vector{C}
\cdot\left(\vector{A}\times\vector{B}\right)=\vector{B}\cdot\left(\vector{C}
\times\vector{A}\right)$:
\begin{equation}
\begin{aligned}\vector{A}\cdot\left(\vector{B}\times\vector{C}\right) & 
=\epsilon_{ijk}A_{i}B_{j}C_{k} & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. 
\ref{EqScalarTripleProduct})}\\
 & =\epsilon_{kij}A_{i}B_{j}C_{k} &  & \text{(Eq. \ref{EqEpsilonCycle})}\\
 & =\epsilon_{kij}C_{k}A_{i}B_{j} &  & \text{(commutativity)}\\
 & =\vector{C}\cdot\left(\vector{A}\times\vector{B}\right) &  & \text{(Eq. 
\ref{EqScalarTripleProduct})}\\
 & =\epsilon_{jki}A_{i}B_{j}C_{k} &  & \text{(Eq. \ref{EqEpsilonCycle})}\\
 & =\epsilon_{jki}B_{j}C_{k}A_{i} &  & \text{(commutativity)}\\
 & =\vector{B}\cdot\left(\vector{C}\times\vector{A}\right) &  & \text{(Eq. 
\ref{EqScalarTripleProduct})}
\end{aligned}
\end{equation}
The negative permutations of these identities can be similarly obtained
and proved by changing the order of the vectors in the cross products
which results in a sign change.

 
$\vector{A}\times\left(\vector{B}\times\vector{C}\right)=\vector{B}\left(\mathbf
{A}\cdot\vector{C}\right)-\vector{C}\left(\vector{A}\cdot\vector{B}\right)$:
\begin{equation}
\begin{aligned}\left[\vector{A}\times\left(\vector{B}\times\vector{C}
\right)\right]_{i} & 
=\epsilon_{ijk}A_{j}\left[\vector{B}\times\vector{C}\right]_{k} & \,\,\,\,\,\, & 
\text{(Eq. \ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}A_{j}\epsilon_{klm}B_{l}C_{m} &  & \text{(Eq. 
\ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\epsilon_{klm}A_{j}B_{l}C_{m} &  & \text{(commutativity)}\\
 & =\epsilon_{ijk}\epsilon_{lmk}A_{j}B_{l}C_{m} &  & \text{(Eq. 
\ref{EqEpsilonCycle})}\\
 & =\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)A_{j}B_{l}C_{m} &  
& \text{(Eq. \ref{EqEpsilonDelta})}\\
 & =\delta_{il}\delta_{jm}A_{j}B_{l}C_{m}-\delta_{im}\delta_{jl}A_{j}B_{l}C_{m} 
&  & \text{(distributivity)}\\
 & 
=\left(\delta_{il}B_{l}\right)\left(\delta_{jm}A_{j}C_{m}\right)-\left(\delta_{
im}C_{m}\right)\left(\delta_{jl}A_{j}B_{l}\right) &  & \text{(commutativity and 
grouping)}\\
 & =B_{i}\left(A_{m}C_{m}\right)-C_{i}\left(A_{l}B_{l}\right) &  & \text{(Eq. 
\ref{EqIndexReplace})}\\
 & 
=B_{i}\left(\vector{A}\cdot\vector{C}\right)-C_{i}\left(\vector{A}\cdot\vector{B
}\right) &  & \text{(Eq. \ref{eqDotProduct})}\\
 & 
=\left[\vector{B}\left(\vector{A}\cdot\vector{C}\right)\right]_{i}-\left[\mathbf
{C}\left(\vector{A}\cdot\vector{B}\right)\right]_{i} &  & \text{(definition of 
index)}\\
 & 
=\left[\vector{B}\left(\vector{A}\cdot\vector{C}\right)-\vector{C}\left(\vector{
A}\cdot\vector{B}\right)\right]_{i} &  & \text{(Eq. \ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.
Other variants of this identity {[}e.g. 
$\left(\vector{A}\times\vector{B}\right)\times\vector{C}${]}
can be obtained and proved similarly by changing the order of the
factors in the external cross product with adding a minus sign.

 
$\vector{A}\times\left(\nabla\times\vector{B}\right)=\left(\nabla\vector{B}
\right)\cdot\vector{A}-\vector{A}\cdot\nabla\vector{B}$:
\begin{equation}
\begin{aligned}\left[\vector{A}\times\left(\nabla\times\vector{B}\right)\right]_
{i} & =\epsilon_{ijk}A_{j}\left[\nabla\times\vector{B}\right]_{k} & 
\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}A_{j}\epsilon_{klm}\partial_{l}B_{m} &  & \text{(Eq. 
\ref{EqCurl})}\\
 & =\epsilon_{ijk}\epsilon_{klm}A_{j}\partial_{l}B_{m} &  & 
\text{(commutativity)}\\
 & =\epsilon_{ijk}\epsilon_{lmk}A_{j}\partial_{l}B_{m} &  & \text{(Eq. 
\ref{EqEpsilonCycle})}\\
 & 
=\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)A_{j}\partial_{l}B_{m
} &  & \text{(Eq. \ref{EqEpsilonDelta})}\\
 & 
=\delta_{il}\delta_{jm}A_{j}\partial_{l}B_{m}-\delta_{im}\delta_{jl}A_{j}
\partial_{l}B_{m} &  & \text{(distributivity)}\\
 & =A_{m}\partial_{i}B_{m}-A_{l}\partial_{l}B_{i} &  & \text{(Eq. 
\ref{EqIndexReplace})}\\
 & =\left(\partial_{i}B_{m}\right)A_{m}-A_{l}\left(\partial_{l}B_{i}\right) &  & 
\text{(commutativity \& grouping)}\\
 & 
=\left[\left(\nabla\vector{B}\right)\cdot\vector{A}\right]_{i}-\left[\vector{A}
\cdot\nabla\vector{B}\right]_{i} &  & \\
 & 
=\left[\left(\nabla\vector{B}\right)\cdot\vector{A}-\vector{A}\cdot\nabla\mathbf
{B}\right]_{i} &  & \text{(Eq. \ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.

 
$\nabla\times\left(\nabla\times\vector{A}\right)=\nabla\left(\nabla\cdot\vector{
A}\right)-\nabla^{2}\vector{A}$:
\begin{equation}
\begin{aligned}\left[\nabla\times\left(\nabla\times\vector{A}\right)\right]_{i} 
& =\epsilon_{ijk}\partial_{j}\left[\nabla\times\vector{A}\right]_{k} & 
\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\partial_{j}\left(\epsilon_{klm}\partial_{l}A_{m}\right) &  & 
\text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\epsilon_{klm}\partial_{j}\left(\partial_{l}A_{m}\right) &  & 
\text{(\ensuremath{\partial} not acting on \ensuremath{\epsilon})}\\
 & =\epsilon_{ijk}\epsilon_{lmk}\partial_{j}\partial_{l}A_{m} &  & \text{(Eq. 
\ref{EqEpsilonCycle} \& definition of derivative)}\\
 & 
=\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)\partial_{j}\partial_
{l}A_{m} &  & \text{(Eq. \ref{EqEpsilonDelta})}\\
 & 
=\delta_{il}\delta_{jm}\partial_{j}\partial_{l}A_{m}-\delta_{im}\delta_{jl}
\partial_{j}\partial_{l}A_{m} &  & \text{(distributivity)}\\
 & =\partial_{m}\partial_{i}A_{m}-\partial_{l}\partial_{l}A_{i} &  & \text{(Eq. 
\ref{EqIndexReplace})}\\
 & =\partial_{i}\left(\partial_{m}A_{m}\right)-\partial_{ll}A_{i} &  & 
\text{(\ensuremath{\partial}\ shift, grouping \& Eq. \ref{eqLaplacianSymbol})}\\
 & 
=\left[\nabla\left(\nabla\cdot\vector{A}\right)\right]_{i}-\left[\nabla^{2}
\vector{A}\right]_{i} &  & \text{(Eqs. \ref{eqDiv}, \ref{eqGrad} \& 
\ref{eqLaplacian2})}\\
 & 
=\left[\nabla\left(\nabla\cdot\vector{A}\right)-\nabla^{2}\vector{A}\right]_{i} 
&  & \text{(Eqs. \ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.
This identity can also be considered as an instance of the identity
before the last one, observing that in the second term on the right
hand side the Laplacian should precede the vector, and hence no independent
proof is required.

 
$\nabla\left(\vector{A}\cdot\vector{B}\right)=\vector{A}
\times\left(\nabla\times\vector{B}\right)+\vector{B}
\times\left(\nabla\times\vector{A}\right)+\left(\vector{A}
\cdot\nabla\right)\vector{B}+\left(\vector{B}\cdot\nabla\right)\vector{A}$:

We start from the right hand side and end with the left hand 
side{\footnotesize{}
\begin{eqnarray}
\left[\vector{A}\times\left(\nabla\times\vector{B}\right)+\vector{B}
\times\left(\nabla\times\vector{A}\right)+\left(\vector{A}
\cdot\nabla\right)\vector{B}+\left(\vector{B}\cdot\nabla\right)\vector{A}\right]
_{i} & =\nonumber \\
\left[\vector{A}\times\left(\nabla\times\vector{B}\right)\right]_{i}+\left[
\vector{B}\times\left(\nabla\times\vector{A}\right)\right]_{i}+\left[
\left(\vector{A}\cdot\nabla\right)\vector{B}\right]_{i}+\left[\left(\vector{B}
\cdot\nabla\right)\vector{A}\right]_{i} & = & \,\,\text{(Eq. 
\ref{eqIndexDistributive1})}\nonumber \\
\epsilon_{ijk}A_{j}\left[\nabla\times\vector{B}\right]_{k}+\epsilon_{ijk}B_{j}
\left[\nabla\times\vector{A}\right]_{k}+\left(A_{l}\partial_{l}\right)B_{i}
+\left(B_{l}\partial_{l}\right)A_{i} & = & \,\,\text{(Eqs. \ref{EqCrossProduct}, 
\ref{eqDiv} \& indexing)}\nonumber \\
\epsilon_{ijk}A_{j}\left(\epsilon_{klm}\partial_{l}B_{m}\right)+\epsilon_{ijk}B_
{j}\left(\epsilon_{klm}\partial_{l}A_{m}\right)+\left(A_{l}\partial_{l}\right)B_
{i}+\left(B_{l}\partial_{l}\right)A_{i} & = & \,\,\text{(Eq. 
\ref{EqCurl})}\nonumber \\
\epsilon_{ijk}\epsilon_{klm}A_{j}\partial_{l}B_{m}+\epsilon_{ijk}\epsilon_{klm}
B_{j}\partial_{l}A_{m}+\left(A_{l}\partial_{l}\right)B_{i}+\left(B_{l}\partial_{
l}\right)A_{i} & = & \,\,\text{(commutativity)}\nonumber \\
\epsilon_{ijk}\epsilon_{lmk}A_{j}\partial_{l}B_{m}+\epsilon_{ijk}\epsilon_{lmk}
B_{j}\partial_{l}A_{m}+\left(A_{l}\partial_{l}\right)B_{i}+\left(B_{l}\partial_{
l}\right)A_{i} & = & \,\,\text{(Eq. \ref{EqEpsilonCycle})}\nonumber \\
\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)A_{j}\partial_{l}B_{m}
+\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)B_{j}\partial_{l}A_{m
}+\left(A_{l}\partial_{l}\right)B_{i}+\left(B_{l}\partial_{l}\right)A_{i} & = & 
\,\,\text{(Eq. \ref{EqEpsilonDelta})}\\
\hspace*{-1.5cm}\left(\delta_{il}\delta_{jm}A_{j}\partial_{l}B_{m}-\delta_{im}
\delta_{jl}A_{j}\partial_{l}B_{m}\right)+\left(\delta_{il}\delta_{jm}B_{j}
\partial_{l}A_{m}-\delta_{im}\delta_{jl}B_{j}\partial_{l}A_{m}\right)+\left(A_{l
}\partial_{l}\right)B_{i}+\left(B_{l}\partial_{l}\right)A_{i} & = & 
\,\,\text{(distributivity)}\nonumber \\
\delta_{il}\delta_{jm}A_{j}\partial_{l}B_{m}-A_{l}\partial_{l}B_{i}+\delta_{il}
\delta_{jm}B_{j}\partial_{l}A_{m}-B_{l}\partial_{l}A_{i}+\left(A_{l}\partial_{l}
\right)B_{i}+\left(B_{l}\partial_{l}\right)A_{i} & = & \,\,\text{(Eq. 
\ref{EqIndexReplace})}\nonumber \\
\delta_{il}\delta_{jm}A_{j}\partial_{l}B_{m}-\left(A_{l}\partial_{l}\right)B_{i}
+\delta_{il}\delta_{jm}B_{j}\partial_{l}A_{m}-\left(B_{l}\partial_{l}\right)A_{i
}+\left(A_{l}\partial_{l}\right)B_{i}+\left(B_{l}\partial_{l}\right)A_{i} & = & 
\,\,\text{(grouping)}\nonumber \\
\delta_{il}\delta_{jm}A_{j}\partial_{l}B_{m}+\delta_{il}\delta_{jm}B_{j}
\partial_{l}A_{m} & = & \,\,\text{(cancellation)}\nonumber \\
A_{m}\partial_{i}B_{m}+B_{m}\partial_{i}A_{m} & = & \,\,\text{(Eq. 
\ref{EqIndexReplace})}\nonumber \\
\partial_{i}\left(A_{m}B_{m}\right) & = & \,\,\text{(product rule)}\nonumber \\
 & = & 
\left[\nabla\left(\vector{A}\cdot\vector{B}\right)\right]_{i}\,\text{(Eqs. 
\ref{eqGrad} \& \ref{eqDiv})}\nonumber
\end{eqnarray}
}Because $i$ is a free index the identity is proved for all components.

 
$\nabla\cdot\left(\vector{A}\times\vector{B}\right)=\vector{B}
\cdot\left(\nabla\times\vector{A}\right)-\vector{A}
\cdot\left(\nabla\times\vector{B}\right)$:
\begin{equation}
\begin{aligned}\nabla\cdot\left(\vector{A}\times\vector{B}\right) & 
=\partial_{i}\left[\vector{A}\times\vector{B}\right]_{i} & 
\,\,\,\,\,\,\,\,\,\,\,\,\,\, & \text{(Eq. \ref{eqDiv})}\\
 & =\partial_{i}\left(\epsilon_{ijk}A_{j}B_{k}\right) &  & \text{(Eq. 
\ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\partial_{i}\left(A_{j}B_{k}\right) &  & 
\text{(\ensuremath{\partial}\ not acting on \ensuremath{\epsilon})}\\
 & =\epsilon_{ijk}\left(B_{k}\partial_{i}A_{j}+A_{j}\partial_{i}B_{k}\right) &  
& \text{(product rule)}\\
 & =\epsilon_{ijk}B_{k}\partial_{i}A_{j}+\epsilon_{ijk}A_{j}\partial_{i}B_{k} &  
& \text{(distributivity)}\\
 & =\epsilon_{kij}B_{k}\partial_{i}A_{j}-\epsilon_{jik}A_{j}\partial_{i}B_{k} &  
& \text{(Eq. \ref{EqEpsilonCycle})}\\
 & 
=B_{k}\left(\epsilon_{kij}\partial_{i}A_{j}\right)-A_{j}\left(\epsilon_{jik}
\partial_{i}B_{k}\right) &  & \text{(commutativity \& grouping)}\\
 & 
=B_{k}\left[\nabla\times\vector{A}\right]_{k}-A_{j}\left[\nabla\times\vector{B}
\right]_{j} &  & \text{(Eq. \ref{EqCurl})}\\
 & 
=\vector{B}\cdot\left(\nabla\times\vector{A}\right)-\vector{A}
\cdot\left(\nabla\times\vector{B}\right) &  & \text{(Eq. \ref{eqDotProduct})}
\end{aligned}
\end{equation}


 
$\nabla\times\left(\vector{A}\times\vector{B}\right)=\left(\vector{B}
\cdot\nabla\right)\vector{A}+\left(\nabla\cdot\vector{B}\right)\vector{A}
-\left(\nabla\cdot\vector{A}\right)\vector{B}-\left(\vector{A}
\cdot\nabla\right)\vector{B}$:
\begin{equation}
\begin{aligned}\hspace{-0.5cm}\left[\nabla\times\left(\vector{A}\times\vector{B}
\right)\right]_{i} & 
=\epsilon_{ijk}\partial_{j}\left[\vector{A}\times\vector{B}\right]_{k} & 
\,\,\,\,\, & \text{(Eq. \ref{EqCurl})}\\
 & =\epsilon_{ijk}\partial_{j}\left(\epsilon_{klm}A_{l}B_{m}\right) &  & 
\text{(Eq. \ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\epsilon_{klm}\partial_{j}\left(A_{l}B_{m}\right) &  & 
\text{(\ensuremath{\partial}\ not acting on \ensuremath{\epsilon})}\\
 & 
=\epsilon_{ijk}\epsilon_{klm}\left(B_{m}\partial_{j}A_{l}+A_{l}\partial_{j}B_{m}
\right) &  & \text{(product rule)}\\
 & 
=\epsilon_{ijk}\epsilon_{lmk}\left(B_{m}\partial_{j}A_{l}+A_{l}\partial_{j}B_{m}
\right) &  & \text{(Eq. \ref{EqEpsilonCycle})}\\
 & 
=\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)\left(B_{m}\partial_{
j}A_{l}+A_{l}\partial_{j}B_{m}\right) &  & \text{(Eq. \ref{EqEpsilonDelta})}\\
 & 
=\delta_{il}\delta_{jm}B_{m}\partial_{j}A_{l}+\delta_{il}\delta_{jm}A_{l}
\partial_{j}B_{m}-\delta_{im}\delta_{jl}B_{m}\partial_{j}A_{l}-\delta_{im}
\delta_{jl}A_{l}\partial_{j}B_{m} &  & \text{(distributivity)}\\
 & 
=B_{m}\partial_{m}A_{i}+A_{i}\partial_{m}B_{m}-B_{i}\partial_{j}A_{j}-A_{j}
\partial_{j}B_{i} &  & \text{(Eq. \ref{EqIndexReplace})}\\
 & 
=\left(B_{m}\partial_{m}\right)A_{i}+\left(\partial_{m}B_{m}\right)A_{i}
-\left(\partial_{j}A_{j}\right)B_{i}-\left(A_{j}\partial_{j}\right)B_{i} &  & 
\text{(grouping)}\\
 & 
=\left[\left(\vector{B}\cdot\nabla\right)\vector{A}\right]_{i}+\left[
\left(\nabla\cdot\vector{B}\right)\vector{A}\right]_{i}-\left[
\left(\nabla\cdot\vector{A}\right)\vector{B}\right]_{i}-\left[\left(\vector{A}
\cdot\nabla\right)\vector{B}\right]_{i} &  & \text{(Eqs. \ref{eqANabla} \& 
\ref{eqDiv})}\\
 & 
=\left[\left(\vector{B}\cdot\nabla\right)\vector{A}+\left(\nabla\cdot\vector{B}
\right)\vector{A}-\left(\nabla\cdot\vector{A}\right)\vector{B}-\left(\vector{A}
\cdot\nabla\right)\vector{B}\right]_{i} &  & \text{(Eq. 
\ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.

 
$\left(\vector{A}\times\vector{B}\right)\cdot\left(\vector{C}\times\vector{D}
\right)=\begin{vmatrix}\begin{array}{cc}
\vector{A}\cdot\vector{C} & \vector{A}\cdot\vector{D}\\
\vector{B}\cdot\vector{C} & \vector{B}\cdot\vector{D}
\end{array}\end{vmatrix}$:
\begin{equation}
\begin{aligned}\left(\vector{A}\times\vector{B}\right)\cdot\left(\vector{C}
\times\vector{D}\right) & 
=\left[\vector{A}\times\vector{B}\right]_{i}\left[\vector{C}\times\vector{D}
\right]_{i} & \,\,\,\,\,\, & \text{(Eq. \ref{eqDotProduct})}\\
 & =\epsilon_{ijk}A_{j}B_{k}\epsilon_{ilm}C_{l}D_{m} &  & \text{(Eq. 
\ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\epsilon_{ilm}A_{j}B_{k}C_{l}D_{m} &  & 
\text{(commutativity)}\\
 & 
=\left(\delta_{jl}\delta_{km}-\delta_{jm}\delta_{kl}\right)A_{j}B_{k}C_{l}D_{m} 
&  & \text{(Eqs. \ref{EqEpsilonCycle} \& \ref{EqEpsilonDelta})}\\
 & 
=\delta_{jl}\delta_{km}A_{j}B_{k}C_{l}D_{m}-\delta_{jm}\delta_{kl}A_{j}B_{k}C_{l
}D_{m} &  & \text{(distributivity)}\\
 & 
=\left(\delta_{jl}A_{j}C_{l}\right)\left(\delta_{km}B_{k}D_{m}
\right)-\left(\delta_{jm}A_{j}D_{m}\right)\left(\delta_{kl}B_{k}C_{l}\right) &  
& \text{(commutativity \& grouping)}\\
 & 
=\left(A_{l}C_{l}\right)\left(B_{m}D_{m}\right)-\left(A_{m}D_{m}\right)\left(B_{
l}C_{l}\right) &  & \text{(Eq. \ref{EqIndexReplace})}\\
 & 
=\left(\vector{A}\cdot\vector{C}\right)\left(\vector{B}\cdot\vector{D}
\right)-\left(\vector{A}\cdot\vector{D}\right)\left(\vector{B}\cdot\vector{C}
\right) &  & \text{(Eq. \ref{eqDotProduct})}\\
 & =\begin{vmatrix}\begin{array}{cc}
\vector{A}\cdot\vector{C} & \vector{A}\cdot\vector{D}\\
\vector{B}\cdot\vector{C} & \vector{B}\cdot\vector{D}
\end{array}\end{vmatrix} &  & \text{(definition of determinant)}
\end{aligned}
\end{equation}


 
$\left(\vector{A}\times\vector{B}\right)\times\left(\vector{C}\times\vector{D}
\right)=\left[\vector{D}\cdot\left(\vector{A}\times\vector{B}\right)\right]
\vector{C}-\left[\vector{C}\cdot\left(\vector{A}\times\vector{B}\right)\right]
\vector{D}$:
\begin{equation}
\begin{aligned}\left[\left(\vector{A}\times\vector{B}\right)\times\left(\vector{
C}\times\vector{D}\right)\right]_{i} & 
=\epsilon_{ijk}\left[\vector{A}\times\vector{B}\right]_{j}\left[\vector{C}
\times\vector{D}\right]_{k} & \,\,\,\,\, & \text{(Eq. \ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\epsilon_{jmn}A_{m}B_{n}\epsilon_{kpq}C_{p}D_{q} &  & 
\text{(Eq. \ref{EqCrossProduct})}\\
 & =\epsilon_{ijk}\epsilon_{kpq}\epsilon_{jmn}A_{m}B_{n}C_{p}D_{q} &  & 
\text{(commutativity)}\\
 & =\epsilon_{ijk}\epsilon_{pqk}\epsilon_{jmn}A_{m}B_{n}C_{p}D_{q} &  & 
\text{(Eq. \ref{EqEpsilonCycle})}\\
 & 
=\left(\delta_{ip}\delta_{jq}-\delta_{iq}\delta_{jp}\right)\epsilon_{jmn}A_{m}B_
{n}C_{p}D_{q} &  & \text{(Eq. \ref{EqEpsilonDelta})}\\
 & 
=\left(\delta_{ip}\delta_{jq}\epsilon_{jmn}-\delta_{iq}\delta_{jp}\epsilon_{jmn}
\right)A_{m}B_{n}C_{p}D_{q} &  & \text{(distributivity)}\\
 & 
=\left(\delta_{ip}\epsilon_{qmn}-\delta_{iq}\epsilon_{pmn}\right)A_{m}B_{n}C_{p}
D_{q} &  & \text{(Eq. \ref{EqIndexReplace})}\\
 & 
=\delta_{ip}\epsilon_{qmn}A_{m}B_{n}C_{p}D_{q}-\delta_{iq}\epsilon_{pmn}A_{m}B_{
n}C_{p}D_{q} &  & \text{(distributivity)}\\
 & =\epsilon_{qmn}A_{m}B_{n}C_{i}D_{q}-\epsilon_{pmn}A_{m}B_{n}C_{p}D_{i} &  & 
\text{(Eq. \ref{EqIndexReplace})}\\
 & =\epsilon_{qmn}D_{q}A_{m}B_{n}C_{i}-\epsilon_{pmn}C_{p}A_{m}B_{n}D_{i} &  & 
\text{(commutativity)}\\
 & 
=\left(\epsilon_{qmn}D_{q}A_{m}B_{n}\right)C_{i}-\left(\epsilon_{pmn}C_{p}A_{m}
B_{n}\right)D_{i} &  & \text{(grouping)}\\
 & 
=\left[\vector{D}\cdot\left(\vector{A}\times\vector{B}\right)\right]C_{i}-\left[
\vector{C}\cdot\left(\vector{A}\times\vector{B}\right)\right]D_{i} &  & 
\text{(Eq. \ref{EqScalarTripleProduct})}\\
 & 
=\left[\left[\vector{D}\cdot\left(\vector{A}\times\vector{B}\right)\right]
\vector{C}\right]_{i}-\left[\left[\vector{C}\cdot\left(\vector{A}\times\vector{B
}\right)\right]\vector{D}\right]_{i} &  & \text{(definition of index)}\\
 & 
=\left[\left[\vector{D}\cdot\left(\vector{A}\times\vector{B}\right)\right]
\vector{C}-\left[\vector{C}\cdot\left(\vector{A}\times\vector{B}\right)\right]
\vector{D}\right]_{i} &  & \text{(Eq. \ref{eqIndexDistributive1})}
\end{aligned}
\end{equation}
Because $i$ is a free index the identity is proved for all components.


\section{The Inertia Tensor}
Consider masses $m_\alpha$ with positions $\vector{r}_\alpha$, all rotating with angular velocity $\omega$ about $\vector{0}$. So the velocities are $\vector{v}_\alpha = \omega\times \vector{r}_\alpha$. The total angular momentum is
\begin{align*}
  \vector{L} &= \sum_{\alpha} \vector{r}_\alpha \times m_\alpha \vector{v}_\alpha \\
  &= \sum_\alpha m_\alpha \vector{r}_\alpha \times (\omega\times \vector{r}_\alpha)\\
  &= \sum_\alpha m_\alpha( |\vector{r}_\alpha|^2\omega - (\vector{r}_\alpha \cdot \omega)\vector{r}_\alpha).
\end{align*}
by vector identities. In components, we have
\[
  L_i = I_{ij}\omega_j,
\]
where
\begin{df}[Inertia tensor]
  The \negrito{inertia tensor} is defined as
  \[
    I_{ij} = \sum_\alpha m_\alpha [|\vector{r}_\alpha|^2 \delta_{ij} - (\vector{r}_\alpha)_i (\vector{r}_\alpha)_j].
  \]
\end{df}
For a rigid body occupying volume $V$ with mass density $\rho(\vector{r})$, we replace the sum with an integral to obtain
\[
  I_{ij} = \dint_V \rho (\vector{r})(x_kx_k \delta_{ij} - x_i x_j)\;\d V.
\]
By inspection, $I$ is a symmetric tensor.

\begin{exa}
  Consider a rotating cylinder with uniform density $\rho_0$. The total mass is $2\ell \pi a^2 \rho_0$.
  \begin{center}
    \begin{tikzpicture}
      \draw [->] (0, 0) -- (2, 0) node [right] {$x_1$};
      \draw [->] (0, 0) -- (0, 2) node [above] {$x_3$};
      \draw [->] (0, 0) -- (-1, -2) node [below] {$x_2$};
      \draw (0.7, -1.5) -- (0.7, 1.5) node [pos = 0.7, right] {$2\ell$};
      \draw (-0.7, -1.5) -- (-0.7, 1.5);
      \draw (0, 1.5) circle [x radius=0.7, y radius=0.3];
      \draw [dashed] (0.7, -1.5) arc (0:180:0.7 and 0.3);
      \draw (-0.7, -1.5) arc (180:360:0.7 and 0.3);
      \draw (0, 1.5) -- (0.7, 1.5) node [pos = 0.5, above] {$a$};
    \end{tikzpicture}
  \end{center}
  Use cylindrical polar coordinate:
  \begin{align*}
    x_1 &= r\cos \theta\\
    x_2 &= r\sin \theta\\
    x_3 &= x_3\\
    \d V &= r\;\d r\;\d \theta \;\d x_3
  \end{align*}
  We have
  \begin{align*}
    I_{33} &= \dint_V \rho_0 (x_1^2 + x_2^2)\;\d V\\
    &= \rho_0 \dint_0^a \dint_0^{2\pi} \dint_{-\ell}^\ell r^2 (r\;\d r\;\d \theta \;\d x_2)\\
    &= \rho_0 \cdot 2\pi \cdot 2\ell \left[\dfrac{r^4}{4}\right]_0^a\\
    &= \varepsilon_0 \pi \ell a^4.
  \end{align*}
  Similarly, we have
  \begin{align*}
    I_{11} &= \dint_V \rho_0 (x_2^2 + x_3^2)\;\d V\\
    &= \rho_0 \dint_0^a \dint_0^{2\pi}\dint_{-\ell}^\ell (r^2 \sin^2 \theta + x_3^2) r\;\d r\;\d \theta \;\d x_3\\
    &= \rho_0 \dint_0^a \dint_0^{2\pi} r\left(r^2 \sin^2 \theta\left[x_3\right]_{-\ell}^\ell + \left[\dfrac{x_3^3}{3}\right]^{\ell}_{-\ell}\right)\;\d \theta\;\d r\\
    &= \rho_0 \dint_0^a \dint_0^{2\pi} r\left(r^2 \sin^2 \theta 2\ell + \dfrac{2}{3}\ell^3\right)\;\d \theta\;\d r\\
    &= \rho_0 \left(2\pi a \cdot \dfrac{2}{3}\ell^3 + 2\ell\dint_0^a r^2 \;\d r\dint_0^{2\pi}\sin^2 \theta\right)\\
    &= \rho_0 \pi a^2 \ell\left(\dfrac{a^2}{2} + \dfrac{2}{3}\ell^2\right)
  \end{align*}
  By symmetry, the result for $I_{22}$ is the same.

  How about the off-diagonal elements?
  \begin{align*}
    I_{13} &= -\dint_V \rho_0 x_1 x_3 \;\d V\\
    &= -\rho_0 \dint_0^a \dint_{-\ell}^\ell \dint_0^{2\pi} r^2 \cos \theta x_3 \;\d r\;\d x_3 \;\d \theta\\
    &= 0
  \end{align*}
  Since $\dint_0^{2\pi} \;\d \theta \cos \theta = 0$. Similarly, the other off-diagonal elements are all 0. So the non-zero components are
  \begin{align*}
    I_{33} &= \dfrac{1}{2}Ma^2\\
    I_{11} = I_{22} &= M\left(\dfrac{a^2}{4} + \dfrac{\ell^2}{3}\right)
  \end{align*}
  In the particular case where $\ell = \dfrac{a\sqrt{3}}{2}$, we have $I_{ij} = \dfrac{1}{2}ma^2 \delta_{ij}$. So in this case,
  \[
    \vector{L} = \dfrac{1}{2}Ma^2 \omega
  \]
  for rotation about any axis.
\end{exa}

\begin{exa}[Inertia Tensor of a Cube about the Center of Mass]

The high degree of symmetry here means we only need to do two out of nine possible integrals. 

\begin{align}
I_{xx} = \dint\ dV \rho (y^2 + z^2) \\  
&= \rho \dint_{-b/2}^{b/2} dx \dint_{-b/2}^{b/2} dy \dint_{-b/2}^{b/2} dz (y^2 + z^2) \\  
&= \rho b \dint_{-b/2}^{b/2} dy \left. (zy^2 + \dfrac{1}{3} z^3) \right|_{-b/2}^{b/2} \\  
&= \rho b \dint_{-b/2}^{b/2} dy \left( by^2 + \dfrac{1}{3} \dfrac{b^3}{4} \right) \\  
&= \rho b \left. \left( \dfrac{1}{3} by^3 + \dfrac{1}{12} b^3 y \right) \right|_{-b/2}^{b/2} \\  
&= \rho b \left( \dfrac{1}{12} b^4 + \dfrac{1}{12} b^4 \right) \\  
&= \dfrac{1}{6} \rho b^5 = \dfrac{1}{6} M b^2.
\end{align}

On the other hand, all the off-diagonal moments are zero, for example
$I_{xy} = \dint\ dV \rho (-xy).$

This is an odd function of x and y, and our integration is now symmetric about the origin in all directions, so it vanishes identically. So the inertia tensor of the cube about its center is

\[\overline{I} = \dfrac{1}{6} Mb^2 \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right).\]
\end{exa}




\subsection{The Parallel Axis Theorem}

The Parallel Axis Theorem relates the inertia  tensor about the center of gravity and the inertia tensor about a parallel axis.

For this purpose we consider two coordinate systems: the first $\vector{r} = (x,y,z)$ with origin at the center of mass of an arbitrary object, and the second $\vector{r'} = (x',y',z')$ offset by some distance. We consider that the object is translated from the origin, but not rotated, by some constant vector $\vector{a}$.

In vector form, the coordinates are related as
\[\vector{r'} = \vector{a} + \vector{r}.\]

Note that $\vector{a}$ points towards the center of mass - the direction is important.
\begin{theorem}
If $I_{ij}$ is the inertia tensor calculated in Center of Mass Coordinate, and $J_{ij}$ is the tensor in the translated coordinates, then:
\[J_{ij} = I_{ij} + M(a^2 \delta_{ij} - a_i a_j). \]
\end{theorem}


\begin{exa}[Inertia Tensor of a Cube about a corner]
The CM inertia tensor was
\[ \overline{I} = Mb^2 \left( \begin{array}{ccc} 1/6 & 0 & 0 \\ 0 & 1/6 & 0 \\ 0 & 0 & 1/6 \end{array} \right)\]

If instead we want the tensor about one corner of the cube, the displacement vector is 
\[\vector{a} = (b/2, b/2, b/2),\]
so $a^2 = (3/4) b^2.$ We can construct the difference as a matrix: the off-diagonal components are
\[M \left[\dfrac{3}{4} B^2 - \left(\dfrac{1}{2} b\right)\left(\dfrac{1}{2} b\right) \right] = \dfrac{1}{2} M b^2\]

and off-diagonal,
\[M \left[-\left(\dfrac{1}{2} b\right)\left(\dfrac{1}{2} b\right) \right] = -\dfrac{1}{4} M b^2\]

so the shifted inertia tensor is
\begin{align}
 \overline{J} &= Mb^2 \left( \begin{array}{ccc} 1/6 & 0 & 0 \\ 0 & 1/6 & 0 \\ 0 & 0 & 1/6 \end{array} \right) + Mb^2 \left( \begin{array}{ccc} 1/2 & -1/4 & -1/4 \\ -1/4 & 1/2 & -1/4 \\ -1/4 & -1/4 & 1/2 \end{array} \right) \\  
&= Mb^2 \left( \begin{array}{ccc} 2/3 & -1/4 & -1/4 \\ -1/4 & 2/3 & -1/4 \\ -1/4 & -1/4 & 2/3 \end{array} \right)
\end{align}


\end{exa}




% 
% \subsection{$\star$ Diagonalization of a symmetric second rank tensor}
% Recall that using matrix notation,
% \[
%   T = (T_{ij}),\quad T' = (T_{ij}'),\quad R = (R_{ij}),
% \]
% and the tensor transformation rule $T'_{ij} = R_{ip}R_{jq}T_{pq}$ becomes
% \[
%   T' = RTR^T = RTR^{-1}.
% \]
% If $T$ is symmetric, it can be diagonalized by such an orthogonal transformation. This means that there exists a basis of orthonormal eigenvectors $\vector{e}_1, \vector{e}_2, \vector{e}_3$ for $T$ with real eigenvalues $\lambda_1, \lambda_2, \lambda_3$ respectively. The directions defined by $\vector{e}_1, \vector{e}_2, \vector{e}_3$ are the \emph{principal axes} for $T$, and the tensor is diagonal in Cartesian coordinates along these axes.
% 
% This applies to any symmetric rank-2 tensor. For the special case of the inertia tensor, the eigenvalues are called the \emph{principal moments of inertia}.
% 
% As exemplified in the previous example, we can often guess the correct principal axes for $I_{ij}$ based on the symmetries of the body. With the axes we chose, $I_{ij}$ was found to be diagonal by direct calculation.
% 
\section{Taylor's Theorem}


\subsection{Multi-index Notation}

An $n$-dimensional \negrito{multi-index} is an $n$-{tuple}

\[\alpha = (\alpha_1, \alpha_2,\ldots,\alpha_n)\]

of natural number.
The set of $n$-{tuple} is denoted \(\mathbb{N}^n_0\).

For multi-indices \(\alpha, \beta \in \mathbb{N}^n_0\) and
\(x = (x_1, x_2, \ldots, x_n) \in \bbR^n\) one defines:

\begin{dingautolist}{202}

\item \textbf{Componentwise sum and difference}
\(\alpha \pm \beta= (\alpha_1 \pm \beta_1,\,\alpha_2 \pm \beta_2, \ldots, 
\,\alpha_n \pm \beta_n)\)

\item \textbf{Partial order}
\(\alpha < \beta \quad \Leftrightarrow \quad \alpha_i < \beta_i \quad 
\forall\,i\in\{1,\ldots,n\}\)

\item \textbf{Sum of components}
\(| \alpha | = \alpha_1 + \alpha_2 + \cdots + \alpha_n\)

\item \textbf{Factorial}
\(\alpha ! = \alpha_1! \cdot \alpha_2! \cdots \alpha_n!\)

\item \textbf{Binomial coefficient}
\(\displaystyle\binom{\alpha}{\beta} = 
\displaystyle\binom{\alpha_1}{\beta_1}\displaystyle\binom{\alpha_2}{\beta_2}\cdots\displaystyle\binom{\alpha_n}{\beta_n
} = \dfrac{\alpha!}{\beta!(\alpha-\beta)!}\)

\item \textbf{Multinomial coefficient}
\(\displaystyle\binom{k}{\alpha} = \dfrac{k!}{\alpha_1! \alpha_2! \cdots \alpha_n! } = 
\dfrac{k!}{\alpha!}\)

\item \textbf{Power}
\(x^\alpha = x_1^{\alpha_1} x_2^{\alpha_2} \ldots x_n^{\alpha_n}\).

\item \textbf{Higher-order Derivatives} \[D^\alpha = D_1^{\alpha_1} 
D_2^{\alpha_2} \ldots D_n^{\alpha_n}\]

\end{dingautolist}
where \(k:=|\alpha|\in\mathbb{N}_0\,\!\) and 
\(\partial_i^{\alpha_i}:=\partial^{\alpha_i} / x_i^{\alpha_i}\)

\subsection{Taylor's Theorem for Multivariate
Functions}\

If the function is $k+1$ times
continuously differentiable in the
closed ball $B$, then one can derive an exact
formula for the remainder in terms of order
partial derivatives of $f$ in this
neighborhood. Namely,

\begin{thm}[Taylor Theorem] \index{Taylor Theorem} Let $f: B\subset \bbR^n \to 
\bbR$ $k+1$ times
continuously differentiable in the
closed ball $B$, then:
\begin{align}& f( \vector{x} ) = \dsum_{|\alpha|<q k} \dfrac{D^\alpha 
f(\vector{a})}{\alpha!}
(\vector{x}-\vector{a})^\alpha  + \dsum_{|\beta|=k+1} 
R_\beta(\vector{x})(\vector{x}-\vector{a})^\beta, \\&
R_\beta( \vector{x} ) = \dfrac{|\beta|}{\beta!} \dint_0^1 (1-t)^{|\beta|-1}D^\beta 
f \big(\vector{a}+t( \vector{x}-
\vector{a} )\big) \, dt. \end{align}

\end{thm}

We prove the special case, where $f:\bbR^n \to \bbR$ has continuous partial
derivatives up to the order $k+1$ in some closed ball $B$ with
center $\vector{a}$.

The strategy of the proof is to apply the
one-variable case of Taylor's theorem to the restriction of $f$ to
the line segment adjoining $\vector{x}$ and
$\vector{a}$. Parametrize the line segment between
$\vector{a}$ and $\vector{x}$ by $u(t) = a + t(x - a)$


We apply the one-variable version
of Taylor's theorem to the function  $g(t) = f(u(t))$

\[f(\vector{x})=g(1)=g(0)+\dsum_{j=1}^k\dfrac{1}{j!}g^{(j)}(0)\ +\ \dint_0^1 
\dfrac{(1-t)^k }{k!} g^{(k+1)}(t)\, dt.\]

Applying the chain rule for several variables gives

\begin{align}
g^{(j)}(t)&=\dfrac{d^j}{dt^j}f(u(t)) = \dfrac{d^j}{dt^j} 
f(\vector{a}+t(\vector{x}-\vector{a})) \\
&= \dsum_{|\alpha|=j} \binom{j}{\alpha} 
(D^\alpha f) (\vector{a}+t(\vector{x}-\vector{a})) 
(\vector{x}-\vector{a})^\alpha
\end{align}

where $\binom{j}{\alpha}$ is the
multinomial coefficient. Since
\(\dfrac{1}{j!}\binom{j}{\alpha}=\dfrac{1}{\alpha!}\),
we get

\[f(\mathbf x)= f(\mathbf a)+\dsum_{|\alpha|<q k}\dfrac{1}{\alpha!} (D^\alpha f) 
(\mathbf a)(\mathbf x-\mathbf a)^\alpha+\dsum_{|\alpha|=k+1}\dfrac{k+1}{\alpha!}
(\mathbf x-\mathbf a)^\alpha \dint_0^1 (1-t)^k (D^\alpha f)(\mathbf a+t(\mathbf 
x-\mathbf a))\,dt.\]


\section{Ohm's Law}


Ohm's law  is an empirical law that states that there is a linear relationship between the electric current $j$ flowing through a material and the electric field $E$ applied to this material. 
This law can be written as
\[j = \sigma E \label{ohm}\]
where the constant of proportionality $\sigma$ is known as the conductivity (the conductivity is defined as the inverse of resistivity). 

One important consequence of   equation \ref{ohm} is that  the vectors $j$ and $E$ are necessary parallel. 

This law is true for some materials, but not for all. For example, if the medium is made of alternate layers of a conductor and an insulator, then the current can only flow along the layers, regardless of the direction of the electric field. It is useful therefore to have an alternative to equation  in which $j$ and $E$ do not have to be parallel. 

This can be achieved by introducing the \negrito{conductivity tensor}, $\sigma_{ik}$, which relates $j$ and $E$ through the equation:
\[j_i=\sigma_{ik}E_k\]

We note that as $j$ and $E$ are vectors, it follows from the quotient rule that $\sigma_{ik}$ is a tensor.





\section{Equation of Motion for a Fluid: Navier-Stokes Equation}
%%%seee Mathew pág 157
% The development below gets relatively involved algebraically. Nevertheless, the key concepts are straight-forward. They are summarized here.
% 
%     The starting point of the Navier-Stokes equations is the equilibrium equation.
%     The first key step is to partition the stress in the equations into hydrostatic (pressure) and deviatoric constituents.
%     The second step is to relate the deviatoric stress to viscosity in the fluid.
%     The final step is to impose any special cases of interest, usually incompressibility.

\subsection{Stress Tensor}
 The stress tensor consists of nine components $\sigma_{ij}$ that completely define the state of stress at a point inside a material in the deformed state, placement, or configuration. 
 

\[\sigma=
\left[{\begin{matrix}
\sigma _{11} & \sigma _{12} & \sigma _{13} \\
\sigma _{21} & \sigma _{22} & \sigma _{23} \\
\sigma _{31} & \sigma _{32} & \sigma _{33} \\
\end{matrix}}\right]
\]



The stress tensor can be separated into two components. One component is a \negrito{hydrostatic} or \negrito{dilatational} stress that acts to change the volume of the material only; the other is the \negrito{deviator stress} that acts to change the shape only.


$$\left( \begin{matrix} {{\sigma _{11}}} & {{\sigma _{12}}} & {{\sigma _{31}}} \cr {{\sigma _{12}}} & {{\sigma _{22}}} & {{\sigma _{23}}} \cr {{\sigma _{31}}} & {{\sigma _{23}}} & {{\sigma _{33}}} \cr \end{matrix} \right) = \left( \begin{matrix} {{\sigma _H}} & 0 & 0 \cr 0 & {{\sigma _H}} & 0 \cr 0 & 0 & {{\sigma _H}} \cr \end{matrix}  \right) + \left( \begin{matrix} {{\sigma _{11}} - {\sigma _H}} & {{\sigma _{12}}} & {{\sigma _{31}}} \cr {{\sigma _{12}}} & {{\sigma _{22}} - {\sigma _H}} & {{\sigma _{23}}} \cr {{\sigma _{31}}} & {{\sigma _{23}}} & {{\sigma _{33}} - {\sigma _H}} \cr \end{matrix}   \right)$$

% where the hydrostatic stress is given by \({\sigma _H}\) = \({1 \over 3}\)\(\left( {{\sigma _1} + {\sigma _2} + {\sigma _3}} \right)\).





\subsection{Derivation of the Navier-Stokes Equations}

The Navier-Stokes equations can be derived from the conservation and continuity equations and some  properties of fluids. In order to derive the equations of fluid motion, we  will first derive the continuity equation, apply the equation to conservation of mass and momentum, and finally combine the conservation equations with a physical understanding of what a fluid is.

The first assumption is that the motion of a fluid are described with the flow velocity of the fluid:

\begin{df}
 The flow velocity $\vector{v}$ of a fluid is a vector field
\[ \vector{v} =\vector{v} (\vector{x} ,t)\]
which gives the velocity of an element of fluid at a position $\vector{x}$ and time $t$ 
\end{df}

\subsubsection*{Material Derivative}



 A normal derivative is the rate of change of of an  property at a point. For instance, the value $\dfrac{dT}{dt}$ could be the rate of change of temperature at a point $(x,y)$. However, a material derivative is the rate of change of an  property on a particle in a velocity field. It incorporates  two things:
\begin{itemize}
    \item Rate of change of the property, $\dfrac{dL}{dt}$
    \item Change in position of of the particle in the velocity field $\vector{v}$
\end{itemize}


Therefore, the material derivative can be defined as

\begin{df}[Material Derivative] Given a function $u(t,x,y,z)$ 
 \[\dfrac{Du}{Dt} = \dfrac{du}{dt} + (\vector{v} \cdot \nabla)u.\]
\end{df}
 



\subsubsection*{Continuity Equation}

An \textbf{intensive property} is a  quantity whose value does not depend on the amount of the substance for which it is measured. For example, the temperature of a system  is the same as the temperature of any part of it. If the system is divided the temperature of each subsystem is identical. The same applies to the density of a homogeneous system; if the system is divided in half, the mass and the volume change in the identical ratio and the density remains unchanged.

The volume will be denoted by  $U$ and its bounding surface area is referred to as $\partial U$. The continuity equation derived can later be applied to mass and momentum.

\paragraph{Reynold's Transport Theorem}

The first basic assumption is the  Reynold's Transport Theorem:


\begin{theorem}[Reynold's Transport Theorem]
 Let \( U\) be a region in $\bbR^n$ with a $C^1$ boundary
\(\partial  U\). Let \(\vector{x}(t)\) be the  positions of
points in the region and let \(\vector{v}(\vector{x},t)\) be the
velocity field in the region. Let \(\vector{n}(\vector{x},t)\) be the
outward unit normal to
the boundary. Let \(L(\vector{x},t)\) be a $C^2$ scalar field. Then
\[\dfrac{d}{dt}\left(\int_{ U} L~\dV\right) = 
     \int_{ U} \frac{\partial L}{\partial t}~\dV + \int_{\partial  U} (\vector{v}\cdot\vector{n})L~\dA ~.\]
\end{theorem}


What we will write in a simplified way as
\begin{equation}
\dfrac{d}{dt}\dint_U L \;\dV = - \dint_{\partial U} L\vector{v} \cdot \vector{n}\;\dA - \dint_U Q\; \dV.
\end{equation}

The left hand side of the equation denotes the rate of change of the property $L$ contained inside the volume $ U$. The right hand side is the sum of two terms:
\begin{itemize}
    \item A flux term, $\dint_{\partial U} L\vector{v} \cdot \vector{n}\;\dA$, which indicates how much of the property $L$ is leaving the volume by flowing over the boundary $\partial U$
    \item A sink term, $\dint_ U Q\; dV$, which describes how much of the property $L$ is leaving the volume due to sinks or sources inside the boundary
\end{itemize}

This equation states that the change in the total amount of a property is due to how much flows out through the volume boundary as well as how much is lost or gained through sources or sinks inside the boundary.


If the intensive property we're dealing with is density, then the equation is simply a statement of conservation of mass: the change in mass is the sum of what leaves the boundary and what appears within it; no mass is left unaccounted for.


% 

% 
% 
% 
% \begin{proof}
% 
% Let \( U _0\) be reference configuration of the region \( U (t)\).
% Let the motion and the deformation gradient be given by
% 
% \[\vector{x} = \vector{\varphi}(\vector{X}, t)~; \qquad\implies\qquad 
%    \vector{F}(\vector{X},t) = \vector{\nabla}_{\circ} \vector{\varphi} ~.\]
% Let \(J(\vector{X},t) = \det[\vector{F}(\vector{X},t)]\). Then,
% integrals in the current and the reference configurations are related by
% 
% \[\int_{ U (t)} L(\vector{x},t)~\dV = 
%       \int_{ U _0} L[\vector{\varphi}(\vector{X},t),t]~J(\vector{X},t)~\dV_0 =
%       \int_{ U _0} \hat{L}(\vector{X},t)~J(\vector{X},t)~\dV_0 ~.\]
% The time derivative of an integral over a volume is defined as
% 
% \[\dfrac{d}{dt}\left( \int_{ U (t)} L(\vector{x},t)~\dV\right) = 
%     \lim_{\Delta t \rightarrow 0} \dfrac{1}{\Delta t}
%      \left(\int_{ U (t + \Delta t)} L(\vector{x},t+\Delta t)~\dV - 
%            \int_{ U (t)} L(\vector{x},t)~\dV\right) ~.\]
% Converting into integrals over the reference configuration, we get
% 
% \[\dfrac{d}{dt}\left( \int_{ U (t)} L(\vector{x},t)~\dV\right) = 
%     \lim_{\Delta t \rightarrow 0} \dfrac{1}{\Delta t}
%      \left(\int_{ U _0} \hat{L}(\vector{X},t+\Delta t)~J(\vector{X},t+\Delta t)~\dV_0 - 
%            \int_{ U _0} \hat{L}(\vector{X},t)~J(\vector{X},t)~\dV_0\right) ~.\]
% Since \( U _0\) is independent of time, we have
% 
% \begin{align}
%   \dfrac{d}{dt}\left( \int_{ U (t)} L(\vector{x},t)~\dV\right) & = 
%     \int_{ U _0} \left[\lim_{\Delta t \rightarrow 0} \dfrac{ 
%            \hat{L}(\vector{X},t+\Delta t)~J(\vector{X},t+\Delta t) - 
%            \hat{L}(\vector{X},t)~J(\vector{X},t)}{\Delta t} \right]~\dV_0 \\
%     & = \int_{ U _0} \frac{\partial }{\partial t}[\hat{L}(\vector{X},t)~J(\vector{X},t)]~\dV_0 \\
%     & = \int_{ U _0} \left(
%           \frac{\partial }{\partial t}[\hat{L}(\vector{X},t)]~J(\vector{X},t)+
%           \hat{L}(\vector{X},t)~\frac{\partial }{\partial t}[J(\vector{X},t)]\right) ~\dV_0 
%   \end{align} Now, the time derivative of \(\det\vector{F}\) 
%   
%   
% \[\frac{\partial J(\vector{X},t)}{\partial t} = \frac{\partial }{\partial t}(\det\vector{F}) = (\det\vector{F})(\vector{\nabla} \cdot \vector{v}) 
%       = J(\vector{X},t)~\vector{\nabla} \cdot \vector{v}(\vector{\varphi}(\vector{X},t),t) 
%       = J(\vector{X},t)~\vector{\nabla} \cdot \vector{v}(\vector{x},t) ~.\]
% Therefore,
% 
% \begin{align}
%   \dfrac{d}{dt}\left( \int_{ U (t)} L(\vector{x},t)~\dV\right) & = 
%      \int_{ U _0} \left(
%           \frac{\partial }{\partial t}[\hat{L}(\vector{X},t)]~J(\vector{X},t)+
%           \hat{L}(\vector{X},t)~J(\vector{X},t)~\vector{\nabla} \cdot \vector{v}(\vector{x},t)\right) ~\dV_0 \\
%      & = 
%      \int_{ U _0} 
%           \left(\frac{\partial }{\partial t}[\hat{L}(\vector{X},t)]+
%           \hat{L}(\vector{X},t)~\vector{\nabla} \cdot \vector{v}(\vector{x},t)\right)~J(\vector{X},t) ~\dV_0  \\
%      & = 
%      \int_{ U (t)} 
%           \left(\dot{L}(\vector{x},t)+
%           L(\vector{x},t)~\vector{\nabla} \cdot \vector{v}(\vector{x},t)\right)~\dV 
%   \end{align} where \(\dot{L}\) is the material time
% derivative of \(L\). Now, the material derivative is given by
% 
% \[\dot{L}(\vector{x},t) = 
%     \frac{\partial L(\vector{x},t)}{\partial t} + [\vector{\nabla} L(\vector{x},t)]\cdot\vector{v}(\vector{x},t) ~.\]
% Therefore,
% 
% \[\dfrac{d}{dt}\left( \int_{ U (t)} L(\vector{x},t)~\dV\right) = 
%      \int_{ U (t)} 
%        \left(
%          \frac{\partial L(\vector{x},t)}{\partial t} + [\vector{\nabla} L(\vector{x},t)]\cdot\vector{v}(\vector{x},t) +
%          L(\vector{x},t)~\vector{\nabla} \cdot \vector{v}(\vector{x},t)\right)~\dV\]
% or,
% 
% \[\dfrac{d}{dt}\left( \int_{ U (t)} L~\dV\right) = 
%      \int_{ U (t)} 
%        \left(
%          \frac{\partial L}{\partial t} + \vector{\nabla} L\cdot\vector{v} +
%          L~\vector{\nabla} \cdot \vector{v}\right)~\dV ~.\]
% Using the identity
% 
% \[\vector{\nabla} \cdot (\vector{v}\otimes\vector{w}) = \vector{v}(\vector{\nabla} \cdot \vector{w}) + \vector{\nabla}\vector{v}\cdot\vector{w}\]
% we then have
% 
% \[\dfrac{d}{dt}\left( \int_{ U (t)} L~\dV\right) = 
%      \int_{ U (t)} 
%        \left(\frac{\partial L}{\partial t} + \vector{\nabla} \cdot (L\otimes\vector{v})\right)~\dV ~.\]
% Using the divergence theorem and the identity
% \((\vector{a}\otimes\vector{b})\cdot\vector{n} = (\vector{b}\cdot\vector{n})\vector{a}\)
% we have
% 
% \[{
%   \dfrac{d}{dt}\left( \int_{ U (t)} L~\dV\right) = 
%      \int_{ U (t)}\frac{\partial L}{\partial t}~\dV + 
%      \int_{\partial  U (t)}(L\otimes\vector{v})\cdot\vector{n}~\dA
%      = \int_{ U (t)}\frac{\partial L}{\partial t}~\dV + 
%      \int_{\partial  U (t)}(\vector{v}\cdot\vector{n})L~\dA ~.
%   }\]
% \end{proof}



\paragraph{Divergence Theorem}

The Divergence Theorem allows the flux term of the above equation to be expressed as a volume integral. By the Divergence Theorem,
\[\dint_{\partial U} L \vector{v}\cdot\vector{n}\;\dA = \dint_{ U} \nabla\cdot(L\vector{v})\;\dV.\]
Therefore, we can now rewrite our previous equation as 
\[\dfrac{d}{dt}\dint_ U L \;\dV = - \dint_{ U}\left[ \nabla\cdot(L\vector{v}) + Q\right]\; dV.\]


Deriving under the integral sign, we find that
\[\dint_ U \dfrac{d}{dt}L \;\dV = - \dint_{ U} \nabla\cdot(L\vector{v}) + Q\; \dV.\]
Equivalently, 
\[\dint_ U \dfrac{d}{dt}L +  \nabla\cdot(L \vector{v}) + Q\;\dV = 0.\]
This relation applies to any volume $ U$; the only way the above equality remains true for any volume $U$ is if the integrand itself is zero. Thus, we arrive at the differential form of the continuity equation
\[\dfrac{dL}{dt} + \nabla\cdot(L\vector{v}) + Q = 0.\]

\subsubsection*{Conservation of Mass}

Applying the continuity equation to density, we obtain
\[\dfrac{d\rho}{dt} + \nabla\cdot(\rho\vector{v}) + Q = 0.\]
This is the  conservation of mass because we are operating with a constant  volume $ U$. With no sources or sinks of mass $(Q=0)$,
\begin{equation}
\dfrac{d\rho}{dt} + \nabla\cdot(\rho\vector{v}) = 0. \label{eq:consofmass}
\end{equation}
The equation \ref{eq:consofmass}  is called   conversation of mass.\\\\


In certain cases it is useful to simplify it further. For an incompressible fluid, the density is constant. Setting the derivative of density equal to zero and dividing through by a constant $\rho$, we obtain the simplest form of the equation
\[\nabla\cdot\vector{v} = 0.\]


\subsubsection*{Conversation of Momentum}

We start with 
\[\vector{F} = m\vector{ a}.\]
Allowing for the body force $\vector{F}=\vector{ a}$ and substituting density for mass, we get a similar equation
\[\vector{b} = \rho \dfrac{d}{dt} \vector{v}(x,y,z,t).\]


Applying the chain rule to the derivative of velocity, we get
\[\vector{b} = \rho\left(\dfrac{\partial \vector{v}}{\partial t} + \dfrac{\partial \vector{v}}{\partial x}\dfrac{\partial x}{\partial t}  + \dfrac{\partial \vector{v}}{\partial y}\dfrac{\partial y}{\partial t} + \dfrac{\partial \vector{v}}{\partial z}\dfrac{\partial z}{\partial t}\right).\]
Equivalently,
\[\vector{b} = \rho\left(\dfrac{\partial \vector{v}}{\partial t} + \vector{v} \cdot \nabla \vector{v}\right).\]
Substituting the value in parentheses for the definition of a material derivative, we obtain 
\begin{equation}
 \rho\dfrac{D\vector{v}}{Dt} = \vector{b}. \label{eq:mom}
\end{equation}

\subsubsection*{Equations of Motion}

The conservation equations derived above, in addition to a few assumptions about the forces and the behaviour of fluids, lead to the equations of motion for fluids. 

We assume that the body force on the fluid parcels is due to two components, fluid stresses and other, external forces. 
\begin{equation}
 \vector{b} = \nabla \cdot \sigma + \vector{f}. \label{eq:mot}
\end{equation}
Here, $\sigma$ is the stress tensor, and $\vector{f}$ represents external forces. Intuitively, the fluid stress is represented as the divergence of the stress tensor because the divergence is the extent to which the tensor acts like a sink or source; in other words, the divergence of the tensor results in a momentum source or sink, also known as a force. 
For many applications  $\vector{f}$ is the gravity force, but for now we will leave the equation in its most general form. 



\subsubsection*{General Form of the Navier-Stokes Equation}

We divide the  stress tensor $\sigma$  into  the hydrostatic and deviator part. Denoting the stress deviator tensor as $T$, we can make the substitution
\begin{equation}
 \sigma = -pI + T.
\end{equation}


Substituting this into the previous equation, we arrive at the most general form of the Navier-Stokes equation:
\begin{equation}
 \rho\dfrac{D\vector{v}}{Dt} = -\nabla p + \nabla \cdot T + \vector{f}.
\end{equation}