图 5.9 求解 (5.25) 的三对角或 Thomas 算法。

function x = tridiag(a,b,c,d,n)
%Solves a tridiagonal system
% a=diag(1 to n), b=upper diag(1 to n-1), c=lower diag(2 to n), d=RHS, %all vectors of dimension n ; note c(1) and b(n)
are not used
% elimination stage
for i=1:n-1
    b(i)=b(i)/a(i);d(i)=d(i)/a(i);a(i)=1;
    a(i+1)=a(i+1)-c(i+1)*b(i);
    d(i+1)=d(i+1)-c(i+1)*d(i);c(i+1)=0;
end
% back substitution
x=zeros(n,1);
d(n)=d(n)/a(n);a(n)=1;
x(n)=d(n);
for j=n-1:-1:1
    x(j)=d(j)-b(j)*x(j+1);
end
%[a,b,c,d] %remove comment at the beginning of the line to
see final a,b,c,d
end

其中

$$b'_1 = \frac{b_1}{a_1}, \quad d'_1 = \frac{d_1}{a_1}, \quad a'_2 = a_2 - c_2 b'_1 \quad \text{and} \quad d'_2 = d_2 - c_2 d'_1$$

接下来我们消去 $x_2$ :

$$\begin{aligned} x_1 + b'_1 x_2 &= d'_1 \ x_2 + b''_2 x_3 &= d''_2 \ a''_3 x_3 + b_3 x_4 &= d''_3 \ c_4 x_3 + a_4 x_4 + b_4 x_5 &= d_4 \end{aligned}$$

依此类推

其中

$$b''_2 = \frac{b_2}{a'_2}, \quad d''_2 = \frac{d'_2}{a'_2}, \quad a''_3 = a_3 - c_3 b''_2 \quad \text{and} \quad d''_3 = d_3 - c_3 d''_2$$

我们可以继续消去所有变量，直到第 $n$ 个。这样我们就把问题转化为上三角形式，可以用图 5.8 中的步骤求解。图 5.9 展示了一个用于求解 (5.25) 的 MATLAB 函数，称为三对角或 Thomas 算法。该算法的编写方式使得每个带撇号的值在计算出来后即替换前一个值。类似地，带双撇号的值替换带撇号的值。这称为覆盖，可减少在计算机上实现该算法所需的存储量。但应当注意，该算法是为清晰而编写的，并非为了最小存储或最高效率。该算法应用非常广泛；它速度极快且所需存储量极少。同样，根据图 5.9 自己编写代码可以大大加深对该方法的理解。