三对角矩阵是指只有系数矩阵的对角线和紧邻对角线的两排是非零矩阵,其他全是零元素的矩阵.如下:
$$\begin{bmatrix} b_0 & c_0 & & & 0 \\
a_1 & b_1 & c_1 \\
& a_2 & b_2 \\
& & \ddots & \ddots & c_{n-2} \\
0 & & & a_{n-1} & b_{n-1} \\ \end{bmatrix} \cdot
\begin{bmatrix} u_0 \\ u_1 \\ \vdots \\ u_{n-2} \\ u_{n-1} \\ \end{bmatrix} =
\begin{bmatrix} r_0 \\ r_1 \\ \vdots \\ r_{n-2} \\ r_{n-1} \\ \end{bmatrix} $$
一个很直观的想法是，因为第0行和最后一行都只有两个元素，如果我们从第0行开始不断的用上一行消掉下一行的a元素，最后一行被会被消掉\(a_{n-1}\)。然后可以直接得到\(x_{n-1}\)。然后不断地从后往前回代，即可求得所有的未知数。这实际上是高斯消元法的一个特例，但是只需要O(N)次操作。给定的算法如下：
$$ c_i^{′} = \begin{cases}
\frac {c_i}{b_i} & i = 0 \\
\frac {c_i}{b_i – a_i c^{′}_{i-1}} & i = 1, 2, \cdots, n-2 \\
\end{cases} $$ $$ r_i^{′} = \begin{cases} \frac{r_i}{b_i} & i = 0 \\ \frac{r_i – a_i r_{i-1}^{′}}{b_i – a_i c_{i-1}^{′}} & i = 1, 2, \cdots, n-1 \\ \end{cases} $$
然后通过回代，可以得到最终结果。
$$\begin{cases}\begin{aligned}x_{n-1} & = r_{n-1}^{′} \\
x_i & = r_{i}^{′} – r_{i}^{′}x_{i+1} \quad i = n-2, n-3, \cdots, 0\end{aligned}\end{cases}$$
上面的公式里实际上涉及到两个动作，假设现在第\(i-1\)行已经消去了\(a_{i-1}\)元素，然后进行的就是以本行的\(b_{i-1}\)元素为基的归一化，然后用本行乘一个系数去消除下一行的\(a_i\)元素。因为\(b_{i-1}\)已经归一，这个系数理所当然的就是下一行的\(a_i\)元素。然后生成的新的\(b_i\)元素就是\(b_i^n = b_i – a_i c_{i-1}^{′}\)，接着对第i行以\(b_i\)为元进行归一化，生成新的\(c_i\)元素（在做行相减操作的时候，\(c_i\)上方的元素是0，所以\(c_i\)的值没有发生变化。）很明显新的\(c_i\)元素的值是：\(c_i^{′} = \frac {c_i}{b_i – a_i c^{′}_{i-1}}\)。在对系数矩阵进行操作的同时，\(r\)列同样要做相应的行相减，乘归一化系数的操作，实际上求\(r_i^{′}\)的公式才是进行变换的完全版。\(c\)只是因为其上的元素为零，省去了相减的操作。经过这样的变换后，所有的\(a\)元素都被消减为0，所有的\(b\)元素都变为1，所以实际上只需要记录\(c^{′}\)元素列和\(r{′}\)列即可。系数矩阵由每行三个有效元素变成每行两个有效元素，最后一行只有一个有效元素。然后我们就可以从后往前推出所有的未知数。
下面给出相应的代码:
[code title=”tridag.h” language=”cpp”]void tridag(VecDoub_I &a, VecDoub_I &b, VecDoub_I &c, VecDoub_I &r, VecDoub_O &u)
{
Int j,n=a.size();
Doub bet;
VecDoub gam(n);
if (b[0] == 0.0) throw("Error 1 in tridag");
u[0]=r[0]/(bet=b[0]);
for (j=1;j<n;j++) {
gam[j]=c[j-1]/bet;
bet=b[j]-a[j]*gam[j];
if (bet == 0.0) throw("Error 2 in tridag");
u[j]=(r[j]-a[j]*u[j-1])/bet;
}
for (j=(n-2);j>=0;j–)
u[j] -= gam[j+1]*u[j+1];
}[/code]
很明显，我们需要要求\(b_0\)和其他的\(b_i – a_i c^{′}_{i-1}\)不能为0。不得不吐槽，这代码写的真的不适合一个本来要去学算法的人看，不过效率确实应该没问题。
下面复述原书中的话：在tridag中没有选主元。因此，矩阵即使非奇异，tridag也可能失败。但在实践中，不必为此担心。那些能产生三对角线性系统的问题通常都有附加的属性，能保证算法tridag成功。当然有可能遇到因为没有选取主要而导致的数值不稳定的情况，但是实践中，这种不稳定几乎从来没有遇到过。
Visits: 1569