这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要复习相关的数学知识。
Vector
Definition
对物理系统的模拟离不开线性代数的相关知识，因此我们从最基础的向量(vector)开始复习。在图形学中向量一般默认是指三维向量，它表示空间点\(\mathbf{p}\)的坐标：
\[\mathbf{p} = \begin{bmatrix} p_x \\ p_y \\ p_z \end{bmatrix} \in \mathbb{R}^3\]在大多数的教材和论文中都会默认我们的坐标系是右手系，但需要注意的是在某些库和引擎中可能会使用左手系(如Unity、DirectX等)。
Addition and Subtraction
向量的基本运算是加法和减法，此时我们只需要将对应位置坐标相加减即可：
\[\mathbf{p} \pm \mathbf{q} = \mathbf{q} \pm \mathbf{p} = \begin{bmatrix} p_x \pm q_x \\ p_y \pm q_y \\ p_z \pm q_z \end{bmatrix} \in \mathbb{R}^3\]从几何的角度上看，向量的加减实际上就是平行四边形法则：
我们可以利用向量的加减法来表示点的运动。假设空间中点在0时刻的位置为\(\mathbf{p}\)并以速度\(\mathbf{v}\)进行运动，则\(t\)时刻\(\mathbf{p}\)点在\(t\)时刻的位置为：
\[\mathbf{p} (t) = \mathbf{p} + t \mathbf{v}\]类似地，我们也可以利用向量加减法来表示直线。经过点\(\mathbf{p}\)和点\(\mathbf{q}\)的直线可以表示为：
\[\begin{aligned} \mathbf{p} (t) &= \mathbf{p} + t (\mathbf{q} - \mathbf{p}) \\ &= (1 - t) \mathbf{p} + t \mathbf{q} \end{aligned}\]当\(0 \leq t \leq 1\)时\(\mathbf{p} (t)\)位于线段\(\mathbf{pq}\)上；当\(t \gt 1\)或\(t \lt 0\)时\(\mathbf{p} (t)\)位于线段\(\mathbf{pq}\)外。从这个角度上看，\(\mathbf{p} (t)\)实际上是对\(\mathbf{p}\)和\(\mathbf{q}\)的位置进行了一个线性插值。
Vector Norm
我们可以通过范数(norm)来度量向量的大小，常用的范数包括：
- 2-norm：\(\Vert p \Vert_2 = (p_x^2 + p_y^2 + p_z^2)^\frac{1}{2}\)
- 1-norm：\(\Vert p \Vert_1 = \vert p_x \vert + \vert p_y \vert + \vert p_z \vert\)
- p-norm：\(\Vert p \Vert_p = (\vert p_x \vert^p + \vert p_y \vert^p + \vert p_z \vert^p)^\frac{1}{p}\)
- infinity norm：\(\Vert p \Vert_\infty = \max \{ \vert p_x \vert, \vert p_y \vert, \vert p_z \vert \}\)
通过范数我们就可以度量\(\mathbf{p}\)和\(\mathbf{q}\)之间的距离\(\Vert \mathbf{p} - \mathbf{q} \Vert\)。同时范数为1的向量称为单位向量(unit vector)，对于任意非零的向量我们可以通过规范化(normalization)来将它转换成单位向量：
\[\mathbf{\bar p} = \frac{\mathbf{p}}{\Vert \mathbf{p} \Vert}\]
Dot Product
点乘运算(dot product)或者说内积运算(inner product)将两个向量映射成标量：
\[\begin{aligned} \mathbf{p} \cdot \mathbf{q} &= p_x q_x + p_y q_y + p_z q_z \\ &= \mathbf{p}^T \mathbf{q} \\ &= \Vert \mathbf{p} \Vert \Vert \mathbf{q} \Vert \cos \theta \end{aligned}\]从几何上看向量之间的点乘等价于它们的模长之积乘以它们夹角的余弦：
点乘运算的常用性质如下：
- 交换律：\(\mathbf{p} \cdot \mathbf{q} = \mathbf{q} \cdot \mathbf{p}\)
- 结合律：\(\mathbf{p} \cdot (\mathbf{q} + \mathbf{r}) = \mathbf{q} \cdot \mathbf{p} + \mathbf{q} \cdot \mathbf{r}\)
- 模长：\(\mathbf{p} \cdot \mathbf{p} = \Vert \mathbf{p} \Vert_2^2\)
- 对于非零向量\(\mathbf{p}\)和\(\mathbf{q}\)，\(\mathbf{p} \cdot \mathbf{q} = 0 \iff \mathbf{p} \perp \mathbf{q}\)
通过点乘运算我们可以完成很多相对复杂的任务。假设我们想要计算\(\mathbf{oq}\)在\(\mathbf{v}\)方向上的投影，首先可以利用点乘得到\(\mathbf{oq}\)在\(\mathbf{v}\)方向上投影的长度：
\[\begin{aligned} s &= \Vert \mathbf{q} - \mathbf{o} \Vert \cos \theta \\ &= \frac{\Vert \mathbf{q} - \mathbf{o} \Vert \Vert v \Vert \cos \theta}{\Vert v \Vert} \\ &= \frac{(\mathbf{q} - \mathbf{o}) \cdot \mathbf{v}}{\Vert v \Vert} \\ &= (\mathbf{q} - \mathbf{o}) \cdot \mathbf{\bar v} \end{aligned}\]因此投影坐标为：
\[\mathbf{s} = \mathbf{o} + s \mathbf{\bar v}\]
类似地，假设平面的法向为\(\mathbf{n}\)，我们可以利用点乘来判断\(\mathbf{p}\)点与平面的相对位置关系：
\[s = (\mathbf{p} - \mathbf{c}) \cdot \mathbf{n}\]若\(s \gt 0\)表示\(\mathbf{p}\)点在平面上方；若\(s \lt 0\)表示\(\mathbf{p}\)点在平面下方；若\(s = 0\)则表示\(\mathbf{p}\)点在平面上。
最后一个例子，我们可以利用点乘来判断\(\mathbf{p}\)点在运动过程中是否与圆盘\(\mathbf{c}\)相交。假设\(t\)时刻\(\mathbf{p}\)点到圆心的距离恰好为\(r\)：
\[\Vert \mathbf{p}(t) - \mathbf{c} \Vert^2 = \Vert \mathbf{p} + t \mathbf{v} - \mathbf{c} \Vert^2 = r^2\]整理方程可以得到关于时间\(t\)的二元一次方程：
\[(\mathbf{v} \cdot \mathbf{v}) t^2 + 2 (\mathbf{p} - \mathbf{c}) \cdot \mathbf{v} t + (\mathbf{p} - \mathbf{c}) \cdot (\mathbf{p} - \mathbf{c}) - r^2 = 0\]因此我们就可以通过这个方程根的数目来判断\(\mathbf{p}\)点是否与圆盘相交。
Cross Product
叉乘运算(cross product)也是一种常用的向量乘法，它定义为：
\[\mathbf{r} = \mathbf{p} \times \mathbf{q} = \begin{bmatrix} p_y q_z - p_z q_y \\ p_z q_x - p_x q_z \\ p_x q_y - p_y q_x \\ \end{bmatrix}\]从几何上看，\(\mathbf{r}\)是一个垂直于\(\mathbf{p}\)和\(\mathbf{q}\)的向量，它的模长等于两个向量模长之积乘以它们夹角的正弦：
\[\Vert \mathbf{r} \Vert = \Vert \mathbf{p} \Vert \Vert \mathbf{q} \Vert \sin \theta\]
叉乘运算的常用性质如下：
- 垂直性：\(\mathbf{r} \cdot \mathbf{p} = \mathbf{r} \cdot \mathbf{q} = 0\)
- 反交换性：\(\mathbf{p} \times \mathbf{q} = -\mathbf{q} \times \mathbf{p}\)
- 结合律：\(\mathbf{p} \times (\mathbf{q} + \mathbf{r}) = \mathbf{q} \times \mathbf{p} + \mathbf{q} \times \mathbf{r}\)
- 对于非零向量\(\mathbf{p}\)和\(\mathbf{q}\)，\(\mathbf{p} \times \mathbf{q} = 0 \iff \mathbf{p} \parallel \mathbf{q}\)
叉乘运算的一个重要应用是可以计算三角形的面积和法向。对于点\(\mathbf{x}_0\)，\(\mathbf{x}_1\)和\(\mathbf{x}_2\)构成的三角形我们可以构造出三角形的两条边：
\[\mathbf{x}_{10} = \mathbf{x}_1 - \mathbf{x}_0\] \[\mathbf{x}_{20} = \mathbf{x}_2 - \mathbf{x}_0\]因此三角形的法向为：
\[\mathbf{n} = \frac{\mathbf{x}_{10} \times \mathbf{x}_{20}}{\Vert \mathbf{x}_{10} \times \mathbf{x}_{20} \Vert}\]三角形的面积则可以表示为：
\[\begin{aligned} A &= \frac{1}{2} \Vert \mathbf{x}_{10} \Vert h \\ &= \frac{1}{2} \Vert \mathbf{x}_{10} \Vert \Vert \mathbf{x}_{20} \Vert \sin \theta \\ &= \frac{1}{2} \Vert \mathbf{x}_{10} \times \mathbf{x}_{20} \Vert \end{aligned}\]不过需要注意的是法向的朝向与三角形三个顶点的顺序有关，因此三角形顶点的顺序也称为拓扑顺序。
类似地，我们可以利用叉乘来判断共面\(\mathbf{p}\)点是否在三角形内部：
当\(\mathbf{p}\)点在三角形内部时会将三角形分成3个部分，它们的面积分别为：
\[A_0 = \frac{1}{2} (\mathbf{x}_1 - \mathbf{p}) \times (\mathbf{x}_2 - \mathbf{p}) \cdot \mathbf{n}\] \[A_1 = \frac{1}{2} (\mathbf{x}_2 - \mathbf{p}) \times (\mathbf{x}_0 - \mathbf{p}) \cdot \mathbf{n}\] \[A_2 = \frac{1}{2} (\mathbf{x}_0 - \mathbf{p}) \times (\mathbf{x}_1 - \mathbf{p}) \cdot \mathbf{n}\] \[A_0 + A_1 + A_2 = A\]我们可以使用面积\(A\)进行归一化，得到：
\[b_0 = \frac{A_0}{A}, b_1 = \frac{A_1}{A}, b_2 = \frac{A_2}{A}\] \[b_0 + b_1 + b_2 = 1\]上式中\((b_0, b_1, b_2)\)称为三角形的重心坐标(barycentric coordinates)，\(\mathbf{p}\)点的坐标可以通过重心坐标插值来得到：
\[\mathbf{p} = b_0 \mathbf{x}_0 + b_1 \mathbf{x}_1 + b_2 \mathbf{x}_2\]
类似地，我们可以利用叉乘来计算四面体的相关几何量：
\[\mathbf{x}_{10} = \mathbf{x}_1 - \mathbf{x}_0\] \[\mathbf{x}_{20} = \mathbf{x}_2 - \mathbf{x}_0\] \[\mathbf{x}_{30} = \mathbf{x}_3 - \mathbf{x}_0\] \[A = \frac{1}{2} \Vert \mathbf{x}_{10} \times \mathbf{x}_{20} \Vert\] \[h = \mathbf{x}_{30} \cdot \mathbf{n} = \mathbf{x}_{30} \cdot \frac{\mathbf{x}_{10} \times \mathbf{x}_{20}}{\Vert \mathbf{x}_{10} \times \mathbf{x}_{20} \Vert}\] \[V = \frac{1}{3} h A = \frac{1}{6} \mathbf{x}_{30} \cdot \mathbf{x}_{10} \times \mathbf{x}_{20} = \frac{1}{6} \left| \begin{matrix} \mathbf{x}_1 & \mathbf{x}_2 & \mathbf{x}_3 & \mathbf{x}_0 \\ 1 & 1 & 1 & 1 \end{matrix} \right|\]
不过需要注意的是按上式计算的体积仍然是带符号的，当\(\mathbf{x}_3\)位于底面三角形上方时为正数，而当\(\mathbf{x}_3\)位于三角形下方时为负数。
类似于三角形的重心坐标，我们同样可以利用体积来定义四面体的重心坐标：
最后来讨论空间射线与三角形求交。我们可以利用\(\mathbf{p}(t)\)点与三角形构成四面体的体积来计算它们的交点，当四面体体积为0时有：
\[(\mathbf{p}(t) - \mathbf{x}_0) \cdot \mathbf{x}_{10} \times \mathbf{x}_{20} = (\mathbf{p} - \mathbf{x}_0 + t \mathbf{v}) \cdot \mathbf{x}_{10} \times \mathbf{x}_{20} = 0\]整理得到：
\[t = \frac{(\mathbf{p} - \mathbf{x}_0) \cdot \mathbf{x}_{10} \times \mathbf{x}_{20}}{\mathbf{v} \cdot \mathbf{x}_{10} \times \mathbf{x}_{20}}\]按上式求解得到的\(t\)即为\(\mathbf{p}\)点到达三角形平面的时刻，接下来只需要检查\(\mathbf{p}(t)\)是否在三角形内部即可。
Matrix
Definition
接下来我们对矩阵(matrix)进行复习，我们把列向量并排放到一起就得到了矩阵。
Multiplication
矩阵的乘法就不再赘述了，这里主要复习一下矩阵乘法的性质：
- 矩阵乘法没有交换律：\(\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}\)
- 矩阵乘法的转置：\((\mathbf{A}\mathbf{B})^T \neq \mathbf{B}^T \mathbf{A}^T\)
- 矩阵(向量)乘以单位阵等于自身：\(\mathbf{I}\mathbf{x} \neq \mathbf{x}\)
- 矩阵的逆：\(\mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I}\)，\((\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}\)
- 不是所有的矩阵都存在逆阵
Orthogonality
由正交的单位向量构成的矩阵称为正交阵(orthogonal matrix)：
\[\mathbf{A} = \begin{bmatrix} \mathbf{a_0} & \mathbf{a_1} & \mathbf{a_2} \end{bmatrix}\] \[\mathbf{a_i}^T \mathbf{a_j} = \left\{ \begin{aligned} 0, & & & i \neq j \\ 1, & & & i = j \\ \end{aligned} \right.\]容易验证正交阵的转置与自身相乘即为单位阵，或者说正交阵的逆等于自身的转置：
\[\mathbf{A}^T \mathbf{A} = \mathbf{I}\] \[\mathbf{A}^{-1} = \mathbf{A}^T\]Matrix Transformation
矩阵在图形学中的一个重要应用是用来表示物体位置和形状的变换。比如说物体的旋转就可以用一个正交阵来表示，假设初始状态物体的三个正方向与世界坐标系相同，旋转后的正方向为\(\mathbf{u}\)，\(\mathbf{v}\)，\(\mathbf{w}\)，则旋转矩阵可以利用这三个方向来表示：
类似地，对物体进行伸缩则可以利用一个对角阵来表示，对角阵上每个元素表示在相应轴上伸缩的比例：
Singular Value Decomposition
对任意矩阵\(\mathbf{A}\)，我们可以对它进行奇异值分解(singular value decomposition, SVD)得到对角阵和正交阵：
\[\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T\]其中\(\mathbf{\Sigma}\)是对角阵，其中的每个元素称为\(\mathbf{A}\)的奇异值；\(\mathbf{U}\)和\(\mathbf{V}\)则是正交阵。
我们可以把矩阵\(\mathbf{A}\)看成是某种线性变换，那么SVD的几何意义则是任意的线性变换都可以通过旋转和缩放来进行描述：
与之类似的还有特征值分解(eigenvalue decomposition)，对称矩阵\(\mathbf{A}\)可以分解为：
\[\mathbf{A} = \mathbf{U} \mathbf{D} \mathbf{U}^{-1}\]其中\(\mathbf{D}\)是对角阵，其中的每个元素称为\(\mathbf{A}\)的特征值；\(\mathbf{U}\)则是正交阵，它的每一列称为\(\mathbf{A}\)的特征向量。矩阵\(\mathbf{A}\)的特征值\(d_i\)和对应的特征向量\(\mathbf{u}_i\)满足：
\[\mathbf{A} \mathbf{u}_i = d_i \mathbf{u}_i\]需要说明的是严格来说特征值分解对于非对称的矩阵仍然成立，但此时的特征值可能是复数。
Symmetric Positive Definiteness
我们称矩阵\(\mathbf{A}\)是一个正定对称阵(symmetric positive definiteness)当且仅当对于任意非零向量\(\mathbf{v}\)满足：
\[\mathbf{v}^T \mathbf{A} \mathbf{v} \gt 0\]称矩阵\(\mathbf{A}\)是一个半正定对称阵(symmetric semi-positive definiteness)当且仅当
\[\mathbf{v}^T \mathbf{A} \mathbf{v} \geq 0\]正定对称矩阵的特征值都大于0，这里简单进行一下证明：
\[\mathbf{v}^T \mathbf{A} \mathbf{v} = \mathbf{v}^T \mathbf{U} \mathbf{D} \mathbf{U}^{T} \mathbf{v} = (\mathbf{U}^{T} \mathbf{v})^T \mathbf{D} (\mathbf{U}^{T} \mathbf{v}) = \sum d_i x_i^2 \gt 0\] \[x_i = \mathbf{U}_i^{T} \mathbf{v}\]由于\(\mathbf{v}\)是任意的，\(\mathbf{v}^T \mathbf{A} \mathbf{v} \gt 0\)要求矩阵\(\mathbf{A}\)的每个特征值都是正数。
实际上特征值都大于0是对称正定矩阵的充要条件，即我们可以利用特征值分解来判断矩阵是否是正定矩阵。但大多数情况下对矩阵进行特征值分解会有较高的计算复杂度，因此我们需要一些简化的方法来进行判断。实践中较为常用的方式是使用对角元素来进行判断，我们称矩阵\(\mathbf{A}\)是对角占优矩阵(diagonally dominant matrix)如果它的对角元素满足：
\[\mathbf{A}_{ii} \gt \sum_{i \neq j} \vert \mathbf{A}_{ij} \vert, i = 1, 2, ..., n\]当一个矩阵是对角占优矩阵时它一定是正定的。但需要注意的是存在非对角占优的正定矩阵。
正定矩阵的一个重要性质是当矩阵\(\mathbf{A}\)是正定矩阵时，矩阵
\[\mathbf{B} = \begin{bmatrix} \mathbf{A} & -\mathbf{A} \\ -\mathbf{A} & \mathbf{A} \end{bmatrix}\]一定是半正定的。这里简单证明一下：
\[\begin{aligned} \begin{bmatrix} \mathbf{x}^T & \mathbf{y}^T \end{bmatrix} \mathbf{B} \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix} &= \begin{bmatrix} \mathbf{x}^T & \mathbf{y}^T \end{bmatrix} \begin{bmatrix} \mathbf{A} & -\mathbf{A} \\ -\mathbf{A} & \mathbf{A} \end{bmatrix} \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix} \\ &= (\mathbf{x} - \mathbf{y})^T \mathbf{A} (\mathbf{x} - \mathbf{y}) \\ &\geq 0 \end{aligned}\]Linear Solver
矩阵的一个重要应用是求解线性方程组：
\[\mathbf{A} \mathbf{x} = \mathbf{b}\]假设矩阵\(\mathbf{A}\)是可逆的方阵，我们可以通过矩阵求逆来求解这个线性方程组。但大多数情况下矩阵求逆具有非常高的计算复杂度，因此实际应用中多半会采取其它的方法来代替矩阵求逆，常用的方法包括直接解法和迭代解法。
直接解法是指对矩阵\(\mathbf{A}\)进行分解来求解线性方程组，常用的矩阵分解包括LU分解、Cholesky分解等。LU分解来计算线性方程组的主要流程如下：
在实践中人们发现当矩阵\(\mathbf{A}\)是稀疏矩阵时直接使用LU分解得到的\(\mathbf{L}\)矩阵和\(\mathbf{U}\)矩阵往往不再是稀疏矩阵，但是可以通过对未知变量\(\mathbf{x}\)进行重新排序来保持矩阵的稀疏性，如在matlab的LU分解中就使用了这样的技巧。
另一种常用的方程解法是迭代解法，我们可以通过建立迭代格式来求解线性方程组：
\[\mathbf{x}^{[k+1]} = \mathbf{x}^{[k]} + \alpha \mathbf{M}^{-1} (\mathbf{b} - \mathbf{A} \mathbf{x}^{[k]})\]其中\(\alpha\)为松弛变量，矩阵\(\mathbf{M}\)为迭代矩阵，\(\mathbf{b} - \mathbf{A} \mathbf{x}^{[k]}\)则是残差向量。可以证明当满足一定条件时残差向量会收敛到\(\mathbf{0}\)：
\[\begin{aligned} \mathbf{b} - \mathbf{A} \mathbf{x}^{[k+1]} &= \mathbf{b} - \mathbf{A} \mathbf{x}^{[k]} - \alpha \mathbf{A} \mathbf{M}^{-1} (\mathbf{b} - \mathbf{A} \mathbf{x}^{[k]}) \\ &= (\mathbf{I} - \alpha \mathbf{A} \mathbf{M}^{-1}) (\mathbf{b} - \mathbf{A} \mathbf{x}^{[k]}) \\ &= (\mathbf{I} - \alpha \mathbf{A} \mathbf{M}^{-1})^{k+1} (\mathbf{b} - \mathbf{A} \mathbf{x}^{[0]}) \end{aligned}\]上式说明当\((\mathbf{I} - \alpha \mathbf{A} \mathbf{M}^{-1})\)的谱半径(spectral radius)小于1时，对于任意初值的变量\(\mathbf{x}\)使用迭代法总能收敛。
除此之外，在使用迭代法求解时还需要计算逆阵\(\mathbf{M}^{-1}\)，这就要求\(\mathbf{M}\)的逆阵需要容易计算。一般可取\(\mathbf{M}\)为：
- Jacobi Method: \(\mathbf{M} = \text{diag}(\mathbf{A})\)
- Gauss-Seidel Method: \(\mathbf{M} = \text{lower} (\mathbf{A})\)
和直接法相比，迭代法的优势在于比较便于实现、容易并行而且往往很快就能收敛到近似解。相对地，使用迭代法时需要注意迭代的收敛条件而且可能需要比较长的时间才能得到精确解。
Tensor Calculus
本节课最后复习了向量和矩阵的微积分。
1st-Order Derivatives
对于标量函数\(f(\mathbf{x})\)，定义它的梯度(gradient)为：
\[\nabla f(\mathbf{x}) = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z} \end{bmatrix}\]梯度\(\nabla f(\mathbf{x})\)的几何意义是在\(\mathbf{x}\)处使向量场增长最快的方向：
对于向量函数\(\mathbf{f} (\mathbf{x}) = \begin{bmatrix} f(\mathbf{x}) \\ g(\mathbf{x}) \\ h(\mathbf{x}) \end{bmatrix}\)可以定义它的Jacobian矩阵为：
\[\mathbf{J} (\mathbf{x}) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}} = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \\ \end{bmatrix}\]其它常用的形式包括散度(divergence)和旋度(curl)：
\[\nabla \cdot \mathbf{f} = \frac{\partial f}{\partial x} + \frac{\partial g}{\partial y} + \frac{\partial h}{\partial z}\] \[\nabla \times \mathbf{f} = \begin{bmatrix} \frac{\partial h}{\partial y} - \frac{\partial g}{\partial z} \\ \frac{\partial f}{\partial z} - \frac{\partial h}{\partial x} \\ \frac{\partial g}{\partial x} - \frac{\partial f}{\partial y} \\ \end{bmatrix}\]2nd-Order Derivatives
对于标量函数\(f(\mathbf{x})\)，我们定义它的二阶导数矩阵为：
\[\mathbf{H} = \mathbf{J} (\nabla f(\mathbf{x})) = \begin{bmatrix} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\ \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial y^2} & \frac{\partial^2 f}{\partial y \partial z} \\ \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial z^2} \\ \end{bmatrix}\]上式称为\(f(\mathbf{x})\)的Hessian矩阵。
另一个常用的二阶微分算子是Laplacian算子，它的定义为：
\[\nabla \cdot \nabla f(\mathbf{x}) = \nabla^2 f(\mathbf{x}) = \Delta f(\mathbf{x}) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2}\]Spring System
利用这些微分算子我们可以计算出很多物理量。假设有一个一端固定的弹簧系统，弹簧的初始长度为\(L\)如下图所示：
当弹簧另一端的位置为\(\mathbf{x}\)时，系统的能量为：
\[E (\mathbf{x}) = \frac{k}{2} (\Vert \mathbf{x} \Vert - L)^2\]弹簧的内力可以通过对能量\(E\)求梯度并取反方向得到：
\[\mathbf{f} (\mathbf{x}) = - \nabla E (\mathbf{x}) = -k (\Vert \mathbf{x} \Vert - L) \frac{\mathbf{x}}{\Vert \mathbf{x} \Vert}\]系统的切向刚度(tangent stiffness)则可以通过对能量\(E (\mathbf{x})\)求二阶导来获得：
\[\mathbf{H} (\mathbf{x}) = - \frac{\mathbf{f} (\mathbf{x})}{\partial \mathbf{x}} = k \frac{\mathbf{x} \mathbf{x}^T}{\Vert \mathbf{x} \Vert^2} + k \bigg(1 - \frac{L}{\Vert \mathbf{x} \Vert} \bigg) \bigg(\mathbf{I} - \frac{\mathbf{x} \mathbf{x}^T}{\Vert \mathbf{x} \Vert^2} \bigg)\]进一步我们可以考虑两端都是自由端的弹簧系统，此时弹簧的伸长量为：
\[\mathbf{x}_{01} = \mathbf{x}_0 - \mathbf{x}_1\]
系统的能量可以表示为：
\[E (\mathbf{x}) = \frac{k}{2} (\Vert \mathbf{x}_{01} \Vert - L)^2\]在这种情况下弹簧系统的能量由\(\mathbf{x}_0\)和\(\mathbf{x}_1\)两个位置矢量控制，对它们求梯度可以得到弹簧的内力：
\[\mathbf{f} (\mathbf{x}) = - \nabla E (\mathbf{x}) = \begin{bmatrix} - \nabla_0 E (\mathbf{x}) \\ - \nabla_1 E (\mathbf{x}) \end{bmatrix} = \begin{bmatrix} \mathbf{f}_e \\ -\mathbf{f}_e \end{bmatrix}\] \[\mathbf{f}_e = -k (\Vert \mathbf{x}_{01} \Vert - L) \frac{\mathbf{x}_{01}}{\Vert \mathbf{x}_{01} \Vert}\]同样地，系统的刚度为：
\[\mathbf{H} (\mathbf{x}) = \begin{bmatrix} \frac{\partial^2 E}{\partial \mathbf{x}_0^2} & \frac{\partial^2 E}{\partial \mathbf{x}_0 \partial \mathbf{x}_1} \\ \frac{\partial^2 E}{\partial \mathbf{x}_0 \partial \mathbf{x}_1} & \frac{\partial^2 E}{\partial \mathbf{x}_1^2} \end{bmatrix} = \begin{bmatrix} \mathbf{H}_e & -\mathbf{H}_e \\ -\mathbf{H}_e & \mathbf{H}_e \end{bmatrix}\] \[\mathbf{H}_e = k \frac{\mathbf{x}_{01} \mathbf{x}_{01}^T}{\Vert \mathbf{x}_{01} \Vert^2} + k \bigg(1 - \frac{L}{\Vert \mathbf{x}_{01} \Vert} \bigg) \bigg(\mathbf{I} - \frac{\mathbf{x}_{01} \mathbf{x}_{01}^T}{\Vert \mathbf{x}_{01} \Vert^2} \bigg)\]