10.2 工程学中的矩阵

一、工程学中的矩阵介绍

这一节介绍的是工程学中的对称矩阵 $K$（通常 $K$ 是正定的）。线性代数中研究对称正定矩阵，我们将其写成 $K=A^{T}A$ 或 $K=A^{T}CA$；物理学中，$\frac{1}{2}{\mathbf{u}}^{T}K\mathbf{u}$ 表示的是能量，能量不可能是负数。其中矩阵 $C$，通常是对角矩阵，其元素表示的是正的物理常数，如电导率（conductance）、刚度（stiffness）或扩散率（diffusivity）。
最好的例子来自于机械工程、土木工程和航空工程。$K$ 是刚度矩阵（stiffness matrix），$K^{-1}\mathbf{f}$ 是系统对外力 $\mathbf{f}$ 的响应。
工程中有两种涉及到线性代数方式，直接的和间接的：
直接方式： 物理问题仅涉及到有限个物体。描述位置和速度的物理定律是线性的（物体的质量不能太大速度不能太快），此时的物理定律可用矩阵方程来描述。
间接方式： 物理系统是 “连续的”。虽然物体的个数有限，但是密度、力和速度是 $x$ 或 $x,y$ 或 $x,y,z$ 的函数，物理定律要用微分方程来表示。为了求得高精度解，我们通常用有限差分方程或有限元方程（finite element equation）来逼近原方程。
这两种方式都会涉及矩阵方程和线性代数。如果没有矩阵，那么在现代工程中什么也做不了。
首先介绍的是平衡方程（equilibrium equation）$K\mathbf{u}=\mathbf{f}$. 如果考虑运动，则会变成动态方程 $M\frac{{\text{d}}^2\mathbf{u}}{\text{d}{t}^2}+K\mathbf{u}=\mathbf{f}$，这是需要使用 $K\mathbf{x}=\lambda M\mathbf{x}$ 对应的特征值或有限差分来求解。

二、从微分方程到矩阵方程

微分方程是连续的，我们的例子是方程 $-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=f(x)$. 矩阵方程是离散的，我们的例子是 $K_0\mathbf{u}=\mathbf{f}$. 通过将二阶导数变成二阶差分，可以将其转换成相应的离散方程。我们先从两个端点 $x=0$ 和 $x=1$ 固定的这样的边值条件开始：

$$\begin{array}{c}\text{Fix-fixedboundaryvalueproblem}\\ 两端固定边值问题\end{array}-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=1,且u(0)=0,u(1)=0(\mathrm{10.2.1})$$

这个微分方程是线性，它的一个特解是 $x_p=-\frac{1}{2}{x}^2$（满足 $\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=-1$），然后加上其 “零空间” 中的任意函数，即求齐次方程 $-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=0$ 的解 $x_n(x)$，它就是 $A\mathbf{x}=0$ 零空间向量 $\mathbf{x}$.
零空间的解是 $u_n(x)=C+Dx$（二阶常微分方程的零空间是二维的），通解（complete solution）就是 $u_p+u_n$：

$$-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=1的通解是u(x)=-\frac{1}{2}{x}^2+C+Dx(\mathrm{10.2.2})$$

下面通过两个边界条件来求解 $C$ 的 $D$：令 $x=0$ 和 $x=1$. 在 $x=0$ 处 $u(0)=0$ 得 $C=0$；在 $x=1$ 处 $u(1)=0$ 得 $-\frac{1}{2}+D=0$，所以 $D=\frac{1}{2}$：$$u(x)=-\frac{1}{2}{x}^2+\frac{1}{2}x=\frac{1}{2}(x-x^2)是这个两端固定边值问题的解.(\mathrm{10.2.3})$$

三、用差分代替导数

要将导数替换为矩阵，我们有三种基本的选择 —— 向前差分、向后差分或中心差分。我们先看一阶导数和一阶差分：$$\frac{\text{d}u}{\text{d}x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x}或\frac{u(x)-u(x-\Delta x)}{\Delta x}或\frac{u(x+\Delta x)-u(x-\Delta x)}{2\Delta x}$$在区间 $x=0$ 和 $x=1$ 之间，我们将其平分成 $n+1$ 分，则每个小区间的长度都是 $\Delta x=\frac{1}{n+1}$. 用差分代替导数后，函数 $u(x)$ 在其 $n$ 个内点 $\Delta x,2\Delta x,\cdots ,n\Delta x$ 的处近似值即为矩阵方程 $K\mathbf{u}=\mathbf{f}$ 中的未知数 $u_1,u_2,\cdots ,u_n$：$$\begin{array}{l}求解：\mathbf{u}=(u_1,u_2,\cdots ,u_n)\approx (u(\Delta x),u(2\Delta x),\cdots ,u(n\Delta x))\\ 边值条件u(0)=u(1)=0得到零值u_0=u_{n+1}=0.\end{array}$$将 $-\frac{\text{d}}{\text{d}x}\left(\frac{\text{d}u}{\text{d}x}\right)=1$ 中的导数替换成向前和向后差分：$$\frac{1}{(\Delta x{)}^2}\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\\ 0& 0& -1\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right](\mathrm{10.2.4})$$这是当 $n=3$ 和 $\Delta x=\frac{1}{4}$ 时得到的矩阵方程，这两个一阶差分矩阵互为转置！这个方式可以写成 $A^{T}A\mathbf{u}=(\Delta x{)}^2\mathbf{f}$，其中 $A$ 表示的是向前差分，$A^T$ 是负的向后差分，注意，这里向前差分的含义是 $A\mathbf{u}=(u_1,u_2-u_1,u_3-u_2,-u_3)$ 表示的是在点 $x=0,1,2,3$ 处一阶微分的近似值，这里包含了边界条件 $u_0=u_4=0$，同理可以得到此条件下 $A^T$ 表示的是负的向后差分矩阵；向前和向后差分的区别是索引不同，如果表示的在是点 $x=1,2,3,4$ 处一阶微分的近似值，则 $A$ 为向后差分矩阵，$A^T$ 为负的向前差分矩阵。由矩阵乘积 $A^{T}A$ 我们得到了一个正定的二阶差分矩阵 $K_0$：

$$K_0\mathbf{u}=(\Delta x{)}^2\mathbf{f}\left[\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 2\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]=\frac{1}{16}\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]得\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]=\frac{1}{32}\left[\begin{array}{c}3\\ 4\\ 3\end{array}\right](\mathrm{10.2.5})$$

这个例子奇妙的地方就是解得的数值 $u_1,u_2,u_3$ 恰好是确切解！它们和真解 $u=\frac{1}{2}(x-x^2)$ 在内点 $x=\frac{1}{4},\frac{2}{4},\frac{3}{4}$ 处的值完全一致。Figure 10.3 展示了真解（连续曲线的解析解）和数值解 $u_1,u_2,u_3$（正好在曲线上），这个曲线是个抛物线。

那么该如何解释数值解恰好落在解析解 $u(x)$ 的图像上这个完美的答案呢？在矩阵方程中，$K_0=A^{T}A$ 是 “二阶差分矩阵”，它给出了 $-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}$ 的中心近似。这里二阶导数包含了负号是因为一阶导数是反对称的，二阶导数本身是负定的：$$\frac{\text{d}}{\text{d}x}的“转置”是-\frac{\text{d}}{\text{d}x}，则(-\frac{\text{d}}{\text{d}x})\left(\frac{\text{d}}{\text{d}x}\right)是正定的.$$令 $A$ 是向前差分矩阵，我们也能够从矩阵 $A$ 和 $A^T$ 中看到这种关系，$A$ 的转置是 $A^T$，它是负的向后差分矩阵。这里没有选择中心差分 $u(x+\Delta x)-u(x-\Delta x)$，因为中心差分对于一阶差分是最好的，但是二阶差分 $A^{T}A$ 将会延伸至 $u(x+2\Delta x)$ 到 $u(x-2\Delta x)$ 这五个函数值，这就很不好。
现在就可以解释上述的完美答案了——差分解恰好落在解析曲线 $u(x)=\frac{1}{2}(x-x^2)$ 上。二阶差分系数 $-1,2,-1$ 对于直线 $y=x$ 和抛物线完全能对应上！$$\begin{array}{lllr}y=x& -\frac{{\text{d}}^{2}y}{\text{d}{x}^2}=0& -(x+\Delta x)+2x-(x-\Delta x)=& 0(\Delta x{)}^2\\ y=x^2& -\frac{{\text{d}}^{2}y}{\text{d}{x}^2}=-2& -(x+\Delta x{)}^2+2x^2-(x-\Delta x{)}^2=& -2(\Delta x{)}^2\end{array}$$这种关系持续到 $y=x^3$，因为 $-\frac{{\text{d}}^{2}y}{\text{d}{x}^2}=-6x$ 也可以通过二阶差分来精确求解。但是对于 $y=x^4$ 将不再成立，二阶差分无法精确对应 $-y^{\prime \prime }=-12x^2$，所以数值解 $u_1,u_2,u_3$ 将不会落在 $u(x)$ 的曲线上。

四、固定端、自由端和变系数 c(x)

为了研究两类新问题，需要改变方程形式和一个边值条件：$$\boxed{-\frac{\text{d}}{\text{d}x}((1+x)\frac{\text{d}u}{\text{d}x})=f(x)且有u(0)=0,\frac{\text{d}u}{\text{d}x}(1)=0}(\mathrm{10.2.6})$$端点 $x=1$ 目前是自由的，此端点没有任何支撑。比如，“一个吊杆只固定了顶部”，在自由端 $x=1$ 处不受任何力，因此在 $x=1$ 处的边值条件不再是固定条件 $u(1)=0$，而变成了 $\frac{\text{d}u}{\text{d}x}(1)=0$.
另一个变化的地方是系数 $c(x)=1+x$，可以表示杆的刚度从 $x=0$ 到 $x=1$ 会发生变化，当杆的宽度或材料变化时会出现该种情况。系数 $1+x$ 将会引入一个新矩阵 $C$ 到差分方程中。
由于 $u_4$ 不再固定为 $0$ 了，它现在是一个新的未知数了。此时的 $A$ 是一个 $4\times 4$ 的向后差分矩阵，引入的系数 $c(x)=1+x$ 则变成了一个对角矩阵 $C$ —— 相当于在分点处分别乘上 $1+\Delta x,1+2\Delta x,1+3\Delta x,1+4\Delta x$. 下面是矩阵 $A^T,C$ 和 $A$：$$A^{T}CA=\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1.25& & & \\ & 1.5& & \\ & & 1.75& \\ & & & 2.0\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ -1& 1& 0& 0\\ 0& -1& 1& 0\\ 0& 0& -1& 1\end{array}\right](\mathrm{10.2.7})$$矩阵 $K=A^{T}CA$ 是对称正定矩阵！由 $(A^{T}CA{)}^T=A^T{C}^T{A}^{TT}=A^{T}CA$ 可推出是对称矩阵。正定可以由能量非负得到：$A$ 的列向量组线性无关，所以当 $\mathbf{x}\ne 0$ 时 $A\mathbf{x}\ne 0$.$$对于任意的\mathbf{x}\ne 0，都有A\mathbf{x}\ne 0，所以：能量={\mathbf{x}}^T{A}^{T}CA\mathbf{x}=(A\mathbf{x}{)}^{T}C(A\mathbf{x})>0$$对于这个固定-自由边值问题，计算乘积 $A^{T}A$ 和 $A^{T}CA$ 时，观察乘积 $A^{T}A$ 右下角的元素是如何从 $2$ 变成 $1$ 的。第四个方程含有 $u_4-u_3$，它来自于自由边值条件 $\frac{\text{d}u}{\text{d}x}(1)=0$ 产生的一阶（不是二阶）差分。
注意 $A^{T}CA$ 中的 $c_1,c_2,c_3,c_4$ 来自于式（10.2.7）中的系数 $c(x)=1+x$，前面两端固定边值的问题中的系数 $c$ 就是简单的 $1,1,1,1$. 下面是固定-自由（fixed-free） 边值问题对应的矩阵：$$A^{T}A=\left[\begin{array}{cccc}2& -1& & \\ -1& 2& -1& \\ & -1& 2& -1\\ & & -1& 1\end{array}\right]A^{T}CA=\left[\begin{array}{cccc}{c}_1+c_2& -c_2& & \\ -c_2& c_2+c_3& -c_3& \\ & -c_3& c_3+c_4& -c_4\\ & & -c_4& c_4\end{array}\right](\mathrm{10.2.8})$$

五、两端自由边值条件

假设杆的两端都是自由的，那么现在在 $x=0$ 和 $x=1$ 处两个端点都有 $\frac{\text{d}u}{\text{d}x}=0$，这种情况下杆没有被固定！物理上它是不稳定的 —— 它可以在没有外力的情况下移动。数学上，所有的常值函数如 $u=1$ 都满足这两端自由的边值条件。代数上，对应的矩阵 $A^{T}A$ 和 $A^{T}CA$ 则是不可逆的：$$\begin{array}{l}两端自由\text{Free-free}的例子\\ 未知数u_0,u_1,u_2,\Delta x=0.5\end{array}{A}^{T}A=\left[\begin{array}{ccc}1& -1& 0\\ -1& 2& -1\\ 0& -1& 1\end{array}\right]A^{T}CA=\left[\begin{array}{ccc}{c}_0& -c_0& \\ -c_0& c_0+c_1& -c_1\\ & -c_1& c_1\end{array}\right]$$其中 $A=\left[\begin{array}{ccc}-1& 1& 0\\ 0& -1& 1\end{array}\right]$，这是因为两端自由缺少了条件，所以 $A$ 是向前差分矩阵，缺少了一行。向量 $(1,1,1)$ 在上述矩阵的零空间中，这与连续模型问题 $u(x)=1$ 的解相吻合。但是对于有外力作用的情形，两端自由问题的方程 $A^{T}A\mathbf{u}=\mathbf{f}$ 和 $A^{T}CA\mathbf{u}=\mathbf{f}$ 通常是不可解的。
在解释更多的物理例子之前，先给出下面的六个矩阵。三对角矩阵 $K_0$ 我们以前经常见到，现在我们看到了它的应用。这些矩阵都是对称的，并且前四个矩阵还是正定的：$$\begin{array}{rr}{K}_0=A_0^T{A}_0=\left[\begin{array}{ccc}2& -1& \\ -1& 2& -1\\ & -1& 2\end{array}\right]& A_0^T{C}_0{A}_0=\left[\begin{array}{ccc}{c}_1+c_2& -c_2& \\ -c_2& c_2+c_3& -c_3\\ & -c_3& c_3+c_4\end{array}\right]\\ 两端固定\text{Fixed-fixed}& 包含弹性系数\text{Spring constants included}\\ K_1=A_1^T{A}_1=\left[\begin{array}{ccc}2& -1& \\ -1& 2& -1\\ & -1& 1\end{array}\right]& A_1^T{C}_1{A}_1=\left[\begin{array}{ccc}{c}_1+c_2& -c_2& \\ -c_2& c_2+c_3& -c_3\\ & -c_3& c_3\end{array}\right]\\ 固定-自由\text{Fiexd-free}& 包含弹性系数\text{Spring constants included}\\ K_{\text{singular}}=\left[\begin{array}{ccc}1& -1& \\ -1& 2& -1\\ & -1& 1\end{array}\right]& K_{\text{circular}}=\left[\begin{array}{ccc}2& -1& -1\\ -1& 2& -1\\ -1& -1& 2\end{array}\right]\\ 两端自由\text{Free-free}& 周期性\text{Periodic}u(0)=u(1)\end{array}$$矩阵 $K_0,K_1,K_{\text{singular}}$ 和 $K_{\text{circular}}$ 均包含简单的系数矩阵 $C=I$，这表明所有的 “弹性系数” 都是 $c_i=1$. 而 $A_0^T{C}_0{A}_0$ 和 $A_1^T{C}_1{A}_1$ 含有弹性系数，这些是为了展示弹性系数如何影响矩阵中的元素，但是不改变其正定性。后面会在其它的工程问题中看到这些刚度矩阵（stiffness matrices）.

六、一列弹簧和质点

Figure 10.4 展示了由一列弹簧（a line of spring：忽略质量）连接起来的三个质点（masses），质量分别是 $m_1,m_2,m_3$. 两端固定的情形有四个弹簧，其顶部和底部都固定，这可以导出 $K_0$ 和 $A_0^T{C}_0{A}_0$. 而固定-自由的情形只有三个弹簧，最下面的质点处于自由悬挂状态，这可以导出 $K_1$ 和 $A_1^T{C}_1{A}_1$. 两端自由的问题会导出矩阵 $K_{\text{singular}}$.
我们希望得出质点的位移 $\mathbf{u}$ 和弹簧张力 $\mathbf{y}$ 的方程：$$\begin{array}{l}\mathbf{u}=(u_1,u_2,u_3)表示每个质点的位移\text{movements of masses}（向下为正）\\ \mathbf{y}=(y_1,y_2,y_3,y_4)或(y_1,y_2,y_3)表示弹簧的张力\text{tensions in the springs}\end{array}$$当质点向下运动时，它的位移为正（$u_j>0$），其中 $u_0$ 和 $u_4$ 表示固定端的位移，由于是固定的，所以有 $u_0=u_4=0$. 对于弹簧，当其伸长时，弹簧的张力会将质点向内拉，张力为正（$y_i>0$），压缩时，弹簧的张力会将质点向外推，张力为负（$y_i<0$）. 每个弹簧都满足其自身的胡克定律（Hooke’s Law）$y=ce$：向量形式为：$\mathbf{y}=\mathbf{c}\mathbf{e}$，其中 $\mathbf{y}$ 表示的是弹力（stretching force），$\mathbf{c}$ 是弹性系数（spring constant），$\mathbf{e}$ 是伸缩距离（stretching distance）.
我们的工作就是要将单个弹簧的方程 $y=ce$ 放到整个系统的向量方程 $K\mathbf{u}=\mathbf{f}$ 中。表示外力的向量 $\mathbf{f}$ 来自于重力，重力加速度常数 $g$ 乘上每个质点的质量得到向下的力 $\mathbf{f}=(m_{1}g,m_{2}g,m_{3}g)$.
问题的关键是要找到刚度矩阵（两端固定和固定-自由这两种情形的）。构造矩阵 $K$ 最好是分成三步，不要一步到位。不要一上来就要建立位移 $u_i$ 和力 $f_i$ 之间的关系，将每个向量按下面的顺序分别和其下面的向量建立联系要更好一些：$$\begin{array}{llll}\mathbf{u}& =& (u_1,u_2,\cdots ,u_n)& 表示n个质点的位移\text{Movements}\\ \mathbf{e}& =& (e_1,e_2,\cdots ,e_m)& 表示m个弹簧的伸缩距离\text{Elongations}\\ \mathbf{y}& =& (y_1,y_2,\cdots ,y_m)& 表示m个弹簧的内部力\text{Internal forces}（即张力）\\ \mathbf{f}& =& (f_1,f_2,\cdots ,f_n)& 表示n个质点受到的外力\text{External forces}\end{array}$$应用数学中会用一个框图将 $\mathbf{u},\mathbf{e},\mathbf{y},\mathbf{f}$ 联系起来，于是得到 $A^{T}CA\mathbf{u}=\mathbf{f}$：

后面将给出这两个例子中的 $A,C$ 和 $A^T$. 首先考虑两端固定的情形，然后再考虑上端固定，下端自由的情形。这些矩阵虽然比较简单，但是它们的形式很重要，尤其是 $A$ 和 $A^T$ 会同时出现。
两端固定情形：$\mathbf{e}$ 表示的是伸缩距离 —— 弹簧的伸缩量。若这个系统是平躺在桌面上时，则弹簧没有伸缩；当垂直放置时，重力就会起作用，质点将会分别下移 $u_1,u_2,u_3$，每个弹簧的伸缩距离为 $e_i=u_i-u_{i-1}$，即是两个端点的位移差：$$\begin{array}{l}每个弹簧的伸缩距离\\ \text{Stretching of each spring}\end{array}\begin{array}{llc}第一个弹簧：& e_1=u_1& (顶端固定，所以u_0=0)\\ 第二个弹簧：& e_2=u_2-u_1& \\ 第三个弹簧：& e_3=u_3-u_2& \\ 第四个弹簧：& e_4=-u_3& (底端固定，所以u_4=0)\end{array}$$如果弹簧的两端都移动相同的距离，那么这个弹簧没有发生伸缩，此时 $u_j=u_{j-1}$，则 $e_j=0$. 上面四个方程得到一个 $4\times 3$ 的差分矩阵 $A$，且有 $\mathbf{e}=A\mathbf{u}$：$$伸缩距离\mathbf{e}=A\mathbf{u}即是\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\\ e_4\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\\ 0& 0& -1\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right](\mathrm{10.2.9})$$下面一个方程 $\mathbf{y}=C\mathbf{e}$ 将弹簧的伸缩距离 $\mathbf{e}$ 和张力 $\mathbf{y}$ 联系了起来，这里对每个弹簧使用了胡克定律 $y_i=c_i{e}_i$. 胡克定律是与弹簧材料有关的 “本构定律 constitutive law”，软弹簧的弹性系数 $c$ 较小，因此强度不大的力 $y$ 即可产生比较大的伸缩距离 $e$. 只要弹簧没有因为被过度拉伸而破坏结构，胡克定律对于弹簧的线性关系是准确的。
因为每个弹簧都遵循自己的定律，所以方程 $\mathbf{y}=C\mathbf{e}$ 中的系数矩阵 $C$ 是一个对角矩阵：$$\begin{array}{c}胡克定律\text{Hooke’s Law}\\ \mathbf{y}=C\mathbf{e}\end{array}\begin{array}{l}{y}_1=c_1{e}_1\\ y_2=c_2{e}_2\\ y_3=c_3{e}_3\\ y_4=c_4{e}_4\end{array}即是\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{cccc}{c}_1& & & \\ & c_2& & \\ & & c_3& \\ & & & c_4\end{array}\right]\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\\ e_4\end{array}\right](\mathrm{10.2.10})$$将方程 $\mathbf{e}=A\mathbf{u}$ 和 $\mathbf{y}=C\mathbf{e}$ 结合起来，可以得到弹簧的张力 $\mathbf{y}=CA\mathbf{u}$.
最后是平衡方程，这是应用数学中最基本的定律。 弹簧张力产生的内部力和质点受到的外力达到平衡。每个质点向上的力来自于其上方弹簧的张力 $y_j$，向下的力是它下方弹簧的张力 $y_{j+1}$ 加上自身的重力 $f_j$. 因此 $y_j=y_{j+1}+f_j$ 或 $f_j=y_j-y_{j+1}$：$$\begin{array}{c}受力平衡\text{Force balance}\\ \mathbf{f}=A^T\mathbf{y}\end{array}\begin{array}{l}{f}_1=y_1-y_2\\ f_2=y_2-y_3\\ f_3=y_3-y_4\end{array}即是\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\end{array}\right](\mathrm{10.2.11})$$这个系数矩阵是 $A^T$！受力平衡方程是 $\mathbf{f}=A^T\mathbf{y}$，将联系 $\mathbf{e}-\mathbf{u}$ 的系数矩阵转置后可以得到联系 $\mathbf{f}-\mathbf{y}$ 的系数矩阵。这个就是框图的美妙之处，$A^T$ 和 $A$ 同时出现了。将这上述的三个方程结合起来就得到 $K\mathbf{u}=\mathbf{f}$.$$\left\{\begin{array}{l}\mathbf{e}=A\mathbf{u}\\ \mathbf{y}=C\mathbf{e}\\ \mathbf{f}=A^T\mathbf{y}\end{array}\right\}\begin{array}{l}组合成A^{T}CA\mathbf{u}=\mathbf{f}或K\mathbf{u}=\mathbf{f}\\ K=A^{T}CA是刚度矩阵\text{stiffness matrix}（力学）\\ K=A^{T}CA是传导矩阵\text{conductance matrix}（网络）\end{array}$$有限元方法（Finite element programs）将需要求解的问题分成多个部分，然后再联系各部分的系数矩阵导出 $K=A^{T}CA$. 我们将 $A^T$ 左乘 $CA$ 就可以得到上述四个弹簧（两端固定）问题的矩阵 $K$：$$\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\end{array}\right]\left[\begin{array}{cccc}{c}_1& & & \\ & c_2& & \\ & & c_3& \\ & & & c_4\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\\ 0& 0& -1\end{array}\right]=\left[\begin{array}{ccc}{c}_1+c_2& -c_2& 0\\ -c_2& c_2+c_3& -c_3\\ 0& -c_3& c_3+c_4\end{array}\right]$$如果所有的弹簧都相同且 $c_1=c_2=c_3=c_4=1$，那么 $C=I$. 则刚度矩阵 $K$ 简化为了 $A^{T}A$，它就变成了特殊的 $-1,2,-1$ 三对角矩阵 $K_0$ 了。
注意区分工程中的 $A^{T}A$ 和线性代数中的 $LU$. 对于四个弹簧的情形矩阵 $A$ 是 $4\times 3$ 的，而消元法中的三角形矩阵是方阵。刚度矩阵 $K$ 是由 $A^{T}A$ 得到的，它可以分解成 $LU$. 前一步是应用数学中的内容，而后一步是计算数学中的内容。应用数学中每个 $K$ 都是由长方阵构造的，后面再分解成方阵的乘积。
下面列出 $K=A^{T}CA$ 的一些性质：

1. $K$ 是三对角矩阵，因为质点 $3$ 和质点 $1$ 不相连；
2. $K$ 是对称矩阵，因为 $C$ 是对称矩阵并且 $A^T$ 是由 $A$ 的转置得来的；
3. $K$ 是正定矩阵，因为 $c_i>0$ 且 $A$ 具有线性无关的列向量组。
4. $K^{-1}$ 是稠密矩阵（不是稀疏矩阵），并且它所有的元素都是正数。

性质 4 揭示了 $\mathbf{u}=K^{-1}\mathbf{f}$ 的一个重要事实：如果所有的力都向下作用（$f_j>0$），则所有的位移都是向下的（$u_j>0$）. 注意 “矩阵元素为正” 与 “矩阵正定” 是不同的。$K^{-1}$ 的所有元素都为正，但是 $K$ 不是，它们都是正定矩阵。

【例1】假设所有的 $c_i=c$ 且 $m_j=m$，求位移 $\mathbf{u}$ 和张力 $\mathbf{y}$.
解： 所有的弹簧都相同，并且所有的质点也相同，但是位移、伸缩距离和张力不会相同。由于 $A^{T}CA$ 含有 $c$，所以 $K^{-1}$ 含有 $\frac{1}{c}$：$$位移\mathbf{u}=K^{-1}\mathbf{f}=\frac{1}{4c}\left[\begin{array}{ccc}3& 2& 1\\ 2& 4& 2\\ 1& 2& 3\end{array}\right]\left[\begin{array}{c}mg\\ mg\\ mg\end{array}\right]=\frac{mg}{c}\left[\begin{array}{c}3/2\\ 2\\ 3/2\end{array}\right]$$由上式可知，中间质点的位移 $u_2$ 要比 $u_1$ 和 $u_3$ 要大。单位也是正确的：力 $mg$（单位：$\text{N}$）除以弹性系数 $c$（单位：$\text{N/m}$），得到 $u$（单位：$\text{m}$）. 则$$伸缩长度\mathbf{e}=A\mathbf{u}=\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\\ 0& 0& -1\end{array}\right]\frac{mg}{c}\left[\begin{array}{c}3/2\\ 2\\ 3/2\end{array}\right]=\frac{mg}{c}\left[\begin{array}{c}3/2\\ 1/2\\ -1/2\\ -3/2\end{array}\right]$$$$\boxed{警告：通常不能写成K^{-1}=A^{-1}{C}^{-1}(A^T{)}^{-1}}$$组成 $A^{T}CA$ 的三个矩阵紧密联系在一起，不太容易能分开。通常，$A^T\mathbf{y}=\mathbf{f}$ 有无穷多解，而这四个方程 $A\mathbf{u}=\mathbf{e}$ 含有三个未知数一般情况下无解。但是在整个系统下，$A^{T}CA$ 可以得到这三个方程的解。当且仅当 $m=n$ 时，这些矩阵都是方阵，我们可以逐步求解 $\mathbf{y}=(A^T{)}^{-1}\mathbf{f}$ 到 $\mathbf{e}=C^{-1}\mathbf{y}$ 到 $\mathbf{u}=A^{-1}\mathbf{e}$. 下面就是这种情形。

七、固定-自由端

移去第四个弹簧，则所有矩阵都会变成 $3\times 3$ 的方阵，但是系统建模的思路不变！此时矩阵 $A$ 不再有第四行，同时 $A^T$ 也不再有第四列。新的刚度矩阵 $K_1$ 就是方阵的乘积：$$A_1^T{C}_1{A}_1=\left[\begin{array}{ccc}1& -1& 0\\ 0& 1& -1\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{c}_1& & \\ & c_2& \\ & & c_3\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\end{array}\right]$$$A^T$ 丢去的列和 $A$ 丢去的行都乘了 $c_4$，所以最快得到新矩阵 $A_1^T{C}_1{A}_1$ 的方法就是在旧的矩阵中令 $c_4=0$：$$固定-自由端A_1^T{C}_1{A}_1=\left[\begin{array}{ccc}{c}_1+c_2& -c_2& 0\\ -c_2& c_2+c_3& -c_3\\ 0& -c_3& c_3\end{array}\right](\mathrm{10.2.12})$$【例2】如果 $c_1=c_2=c_3=1$，则 $C=I$，那么此时矩阵 $K_1$ 就是类似于 “$-1,2,1$” 的三对角矩阵，由于弹簧的底部没有固定，所以 K_1 最后一个元素是 $1$ 而不是 $2$. 设所有的 $m_j=m$，弹性系数均为 $c$：$$固定-自由端\mathbf{u}=K_1^{-1}\mathbf{f}=\frac{1}{c}\left[\begin{array}{ccc}1& 1& 1\\ 1& 2& 2\\ 1& 2& 3\end{array}\right]\left[\begin{array}{c}mg\\ mg\\ mg\end{array}\right]=\frac{mg}{c}\left[\begin{array}{c}3\\ 5\\ 6\end{array}\right]$$这些位移要比两端固定的情况更大。${\mathbf{u}}_1$ 的值是 $\frac{3mg}{c}$ 是因为三个质点都将第一个弹簧往下拉；第二个质点的位移是第一个质点的位移 $\frac{3mg}{c}$ 再额外加上下面的两个质点将第二个弹簧拉长的距离 $\frac{2mg}{c}$；第三个质点下降的最多，其位移是 $\frac{3mg}{c}+\frac{2mg}{c}+\frac{mg}{c}$. 弹簧的伸缩距离 $\mathbf{e}=A\mathbf{u}$，可得其包含了数字 $3,2,1$：$$\mathbf{e}=\left[\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& -1& 1\end{array}\right]\frac{mg}{c}\left[\begin{array}{c}3\\ 5\\ 6\end{array}\right]=\frac{mg}{c}\left[\begin{array}{c}3\\ 2\\ 1\end{array}\right]$$

八、两个自由端：K 是奇异的

如果两端都是自由的，我们会遇到麻烦。因为这整个系统可以移动，此时 $A$ 是 $2\times 3$ 的矩阵：$$\begin{array}{l}两端自由\\ e=A\mathbf{u}\end{array}\left[\begin{array}{c}{e}_1\\ e_2\end{array}\right]=\left[\begin{array}{c}{u}_2-u_1\\ u_3-u_2\end{array}\right]=\left[\begin{array}{ccc}-1& 1& 0\\ 0& -1& 1\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right](\mathrm{10.2.13})$$此时 $A\mathbf{u}=0$ 有一个非零解，质点可以在弹簧不伸缩的情况下运动。如，质点可以移动 $\mathbf{u}=(1,1,1)$，而此时的 $\mathbf{e}=(0,0)$：$$A\mathbf{u}=\left[\begin{array}{ccc}-1& 1& 0\\ 0& -1& 1\end{array}\right]\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]，此时没有伸缩(\mathrm{10.2.14})$$$A\mathbf{u}=0$ 会使得 $A^{T}CA\mathbf{u}=0$，矩阵 $A^{T}CA$ 中没有 $c_1$ 和 $c_4$，它只是一个半正定的矩阵，主元是 $c_2$ 和 $c_3$，没有第三个主元了，秩为 $2$：$$\left[\begin{array}{cc}-1& 0\\ 1& -1\\ 0& 1\end{array}\right]\left[\begin{array}{cc}{c}_2& \\ & c_3\end{array}\right]\left[\begin{array}{ccc}-1& 1& 0\\ 0& -1& 1\end{array}\right]=\left[\begin{array}{ccc}{c}_2& -c_2& 0\\ -c_2& c_2+c_3& -c_3\\ 0& -c_3& c_3\end{array}\right](\mathrm{10.2.15})$$这个矩阵的两个特征值为正，但是 $\mathbf{x}=(1,1,1)$ 是特征值 $\lambda =0$ 的一个特征向量。我们只能针对一些特殊的向量 $\mathbf{f}$ 才能够求解 $A^{T}CA=\mathbf{f}$，外力要满足 $f_1+f_2+f_3=0$，否则整个系统（两个端点自由）将会一直向上或向下运动。

九、弹簧圈

两端自由的情形，再加上第三个弹簧，将质点 $3$ 和质点 $1$ 连接起来可以形成一个圈，这就是一个弹簧圈（circle of springs）. 此时的矩阵 $K$ 仍然是不可逆矩阵 —— 刚度矩阵 $K_{\text{circular}}$ 仍然是奇异的：$$A_{\text{circular}}^T{A}_{\text{circular}}=\left[\begin{array}{ccc}1& -1& 0\\ 0& 1& -1\\ -1& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& -1\\ -1& 1& 0\\ 0& -1& 1\end{array}\right]=\left[\begin{array}{ccc}2& -1& -1\\ -1& 2& -1\\ -1& -1& 2\end{array}\right](\mathrm{10.2.16})$$这个矩阵只有两个主元 $2$ 和 $\frac{3}{2}$，特征值是 $3,3$ 和 $0$，行列式为零。零空间仍然包含 $\mathbf{x}=(1,1,1)$，说明所有的质点会一起移动。 位移向量 $(1,1,1)$ 在矩阵 $A_{\text{circular}}$ 和 $K_{\text{circular}}=A^{T}CA$ 的零空间中。

下面是对本节的总结： 用差分方程近似求解微分方程，将线性代数和微积分联系在了一起。如果所取的步长 $\Delta x$ 足够小，那么就可以获得让自己足够满意的解。$$方程-\frac{\text{d}}{\text{d}x}(c(x)\frac{\text{d}u}{\text{d}x})=f(x)满足u(0)=0且u(1)=0或\frac{\text{d}u}{\text{d}x}(1)=0$$将区间 $[0,1]$ 分成 $N$ 分，每分区间长度为 $\Delta x$，用 $A\mathbf{u}$ 代替 $\frac{\text{d}u}{\text{d}x}$，$A^T\mathbf{y}$ 代替 $-\frac{\text{d}y}{\text{d}x}$（其中 $y=c(x)\frac{\text{d}y}{\text{d}x}$）. 则此时 $A$ 和 $A^T$ 中的元素都含有 $\frac{1}{\Delta x}$. 边值条件是 $u_0=0$ 和 $u_N=0$ 或 $y_N=0$. 这三步 $-\frac{\text{d}}{\text{d}x},c(x)$ 和 $\frac{\text{d}}{\text{d}x}$ 对应的矩阵分别是 $A^T,C$ 和 $A$：$$\mathbf{f}=A^T\mathbf{y},\mathbf{y}=C\mathbf{e},\mathbf{e}=A\mathbf{u}得到A^{T}CA=\mathbf{f}$$这是计算科学可工程中的一个基本的例子，实际建模和求解的步骤如下：

1. 用微分方程描述实际问题；
2. 将微分方程离散成差分方程；
3. 求解差分方程（利用边值条件）；
4. 解释解的含义，绘制图形，必要时要重新建模。

数值模拟已经成为了科学领域的第三个分支，另外两个是实验和理论。比如在计算机上设计波音 $777$ 要比风洞实验成本要低得多。