數值分析（二）：C++實現三對角線方程組的追趕法
這次來實現三對角線方程組的追趕法，追趕法的本質還是高斯消元法，而且是沒選主元的高斯消元法，只是因為Ax=b中係數矩陣A非常特殊，所以就可以採用相對特殊的方法來解方程組。同樣，按照常規的步驟，先分析什麼是追趕法，再給出追趕法的數學步驟，最後用C++實現這種演算法。
（一）追趕法的功能和步驟
明確好目的，正所謂磨刀不誤砍柴工，做一件事情事先規劃好，那重要性真的是不言而喻。在一些實際問題中，對角佔優的三對角線方程組很常見，如熱傳導方程，形式如下圖1：
應用追趕法能解這種三對角線方程組還有很嚴苛的條件，那就是對角佔優，主對角線元素的絕對值要最大，大到什麼程度呢？大到絕對值比旁邊兩條次對角線的值的絕對值之和還要大，這樣才能用追趕法來接三對角線方程組，苛刻的條件如下圖2：
所以，追趕法就是用來解這類線性方程組的。
追趕法的步驟：追趕法的本質就是沒有選主元的高斯消元法，當然係數矩陣***A***都這麼特殊了，再用原始的高斯消元法就實在是浪費空間和時間了，特殊的特性就應該用起來，就像做數學證明題時一樣，要把事物本身的性質儘可能都用上。實現追趕法最有趣的不是數學方法，而是程式設計實現，因為之前的係數矩陣***A***是低階的稠密矩陣，所以程式設計採用的資料結構就是二維的矩陣，或者說是二重指標，但是實現追趕法，沒必要用到二重指標，只要一維陣列就夠用了，因為就三對角線上有元素，所以一維陣列的大小是**n+2(n-1)***，即**3n-1***就夠了，然後最重要的就是將係數矩陣***A**的
所以，資料結構確定，且對應關係理清之後，就又是沒選主元的高斯消元法了。數學語言的步驟如下圖4。
（二）程式碼實現
筆者分別將追趕法寫成了函式和類兩種方式，仔細介紹函式double trde(int n, double A, double b)***，函式就是一個處理器，或是黑匣子，或說是類似對映一樣的東西，輸入引數，計算處理後得到我們想要的輸出，這裡引數就是係數矩陣的階數***n***，以及儲存有係數矩陣元素的一維陣列***A***（當然這裡以指標的形式傳輸），以及常數向量***b***，輸出就是得到解向量，同樣是指標形式。程式碼如下：
double* trde(int n, double* A, double* b)
{
int k, j, m; m = 3 * n - 2;
double s;
for (k = 0; k < n - 1; k++)
{
j = 3 * k; s = A[j];
if (fabs(3) + 1.0 == 1.0)
{
cout << "主對角元素是0，方程無解！！" << endl;
return b;
}
A[j + 1] = A[j + 1] / s;//一下這四行是歸一化和消元，對角線下的元素沒必要進行操作
b[k] = b[k] / s;//因為後面的回代沒用到這些資料
A[j + 3] = A[j + 3] - A[j + 2] * A[j + 1];
b[k + 1] = b[k + 1] - A[j + 2] * b[k];
}
s = A[3 * n - 3];
if (fabs(s) + 1.0 == 1.0)
{
cout << "消元后，主對角元素是0，無解！！" << endl;
return b;
}
b[n - 1] = b[n - 1] / s;
for (k = n - 2; k >= 0; k--)//回代求解方程，相比之下很簡單了，累積求和都沒必要的
b[k] = b[k] - b[k + 1] * A[3 * k + 1];
return b;
}
主函式的呼叫：
int main()
{
int n = 5;
double A[13]{ 2,-1,-1,2,-1,-1,2,-1 ,-1,2,-1 ,-1,2 };
double b[5]{ 1,0,0,0,0 };
double *d;
d = trde(n, A, b);
for (int i = 0; i < n; i++)
cout << d[i] << " ";
//這裡直接將解顯示在黑框框裡了
system("pause");
return 0;
}
執行結果，原線性方程組，以及輸出的解，如圖5,6所示。
最後把類的形式給出來，和全選主元高斯消元法的形式是一樣的，input（），輸入txt檔名，然後進行操作，a_trde（），輸出結果儲存在txt檔案裡，自己命名，接下來是完整的追趕法的類的原始碼。
//追趕求解三對角方程組
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
class trde
{
private:
int n;
double *b, *d;
public:
trde(int nn)
{
n = nn;
b = new double[3 * n - 2]; //動態分配記憶體空間，一維陣列，存放係數矩陣A的元素
d = new double[n];//一維陣列存放常數向量b的元素
}
void input(); //從檔案讀入存放三對角矩陣的向量B以及常數向量D
void a_trde(); //執行追趕法
void output(); //輸出結果到檔案並顯示
~trde()
{
delete[] b;
delete[] d;
}
};
void trde::input() //從檔案讀入存放三對角矩陣的向量B以及常數向量D
{
int i;
char str1[20];
cout << "\n輸入檔名: ";
cin >> str1;
ifstream fin(str1);
if (!fin)
{
cout << "\n不能開啟這個檔案 " << str1 << endl; exit(1);
}
for (i = 0; i<3 * n - 2; i++) fin >> b[i]; //讀入三對角矩陣A中的非0元素
for (i = 0; i<n; i++) fin >> d[i]; //讀入常數向量d
fin.close();
}
void trde::a_trde() //執行追趕法，核心的函式
{
int k, j, m;
double s;
m = 3 * n - 2;//第一步，對於k從0~n-2，執行歸一和消元：
for (k = 0; k <= n - 2; k++)
{
j = 3 * k; s = b[j];
if (fabs(s) + 1.0 == 1.0)
{
cout << "\n程式工作失敗！無解. " << endl;
return;
}
b[j + 1] = b[j + 1] / s;//係數矩陣的歸一化
d[k] = d[k] / s;//常數向量的歸一化
b[j + 3] = b[j + 3] - b[j + 2] * b[j + 1];//係數矩陣的消元
d[k + 1] = d[k + 1] - b[j + 2] * d[k];//消元過程常數向量的變化
}//第二步：進行回代，倒序進行求解
s = b[3 * n - 3];
if (fabs(s) + 1.0 == 1.0)
{
cout << "\n程式工作失敗！無解. " << endl;
return;
}
d[n - 1] = d[n - 1] / s;
for (k = n - 2; k >= 0; k--)
d[k] = d[k] - b[3 * k + 1] * d[k + 1];
}
void trde::output()//結果輸出，自己命名txt檔案，解向量儲存在檔案裡
{
int i;
char str[20];
cout << "解儲存在txt檔案裡，輸入檔名！" << endl;
cin >> str;
ofstream fout(str);
if (!fout)
{
cout << "不能開啟這個檔案！！" << str << endl;
exit(1);
}
for (i = 0; i < n; i++)
{
fout << d[i] << " ";
cout << d[i] << " ";
}
fout << endl; cout << endl;
fout.close();
}
//接下來是主函式的呼叫
/*
int main() //主函式
{
trde c(5);
c.input(); //從檔案讀入存放三對角矩陣的向量B以及常數向量D
c.a_trde(); //執行追趕法
c.output(); //輸出結果到檔案並顯示
system("pause");
return 0;
}
*/
上週六一看數值分析老師的進度，嚇一跳，原來自己已經落後這麼多了（捂臉哭），所以又學習程式設計實現解三對角線方程組的追趕法，再補充一篇，自己本沒有時間觀念，但是也算體會到了動漫團隊持續一週更新一集動漫的艱辛，若是讓我每一週都學習並用C++實現一個演算法，那我肯定堅持不住了。佩服那些堅持周更的作者們，真是不容易，正好前陣子看了鳴人的原畫師作者黃成希的一席演講，真的很不容易啊，那樣的勞動強度，為我們展示瞭如此精彩的動漫熱血視訊，人家用十年時間成為了最優秀的原畫師，想想我十年能夠成為最優秀的某某某，還真是羞愧難當啊，只能去追求美好生活了，哈哈~