考察线性方程组：
$$
Ax=y
$$
其中系数矩阵A为如下形式：

1. 常对角矩阵(Toeplitz矩阵)
   $$
   A = \begin{bmatrix}
   a_{0} & a_{-1} & a_{-2} & \ldots & \ldots & a_{-n+1} \
   a_{1} & a_{0} & a_{-1} & \ddots & & \vdots \
   a_{2} & a_{1} & \ddots & \ddots & \ddots & \vdots \
   \vdots & \ddots & \ddots & \ddots & a_{-1} & a_{-2} \
   \vdots & & \ddots & a_{1} & a_{0} & a_{-1} \
   a_{n-1} & \ldots & \ldots & a_{2} & a_{1} & a_{0}
   \end{bmatrix}
   $$
2. 循环矩阵
   $$
   A = \begin{bmatrix}
   a_{0} & a_{n-1} & a_{n-2} & \ldots & \ldots & a_{1} \
   a_{1} & a_{0} & a_{n-1} & a_{n-2} & & \vdots \
   a_{2} & a_{1} & a_{0} & a_{n-1} & \ddots & \vdots \
   \vdots & \ddots & \ddots & \ddots & a_{1} & a_{n-2} \
   \vdots & & \ddots & a_{1} & a_{0} & a_{n-1} \
   a_{n-1} & a_{n-2} & a_{n-3} & \ldots & \ldots & a_{0}
   \end{bmatrix}
   $$
   为特殊的Toeplitz矩阵。
3. 三对角矩阵
   $$
   A = \begin{bmatrix}
   b_{1} & c_{1} & & & \
   a_{2} & b_{2} & c_{2} & & \
   & a_{3} & b_{3} & \ddots & \
   & & \ddots & \ddots & c_{n-1}\
   & & & a_{n} & b_{n} \
   \end{bmatrix}
   $$
4. 循环Toeplitz三对角线性方程组
   $$
   A = \begin{bmatrix}
   b_{1} & c_{1} & & & a_{1}\
   a_{2} & b_{2} & c_{2} & & \
   & a_{3} & b_{3} & \ddots & \
   & & \ddots & \ddots & c_{n-1}\
   c_{n} & & & a_{n} & b_{n} \
   \end{bmatrix}
   $$
5. 五对角矩阵
   $$
   A = \begin{bmatrix}
   c_{1} & d_{1} & e_{1} & & & & \
   b_{2} & c_{2} & d_{2} & e_{2} & & & \
   a_{3} & b_{3} & c_{3} & d_{3} & e_{3} & & \
   & a_{4} & b_{4} & c_{4} & d_{4} & e_{4} & \
   & & \ddots & \ddots & \ddots & \ddots & \ddots \
   & & & a_{n-2} & b_{n-2} & c_{n-2} & d_{n-2} & e_{n-2} \
   & & & & a_{n-1} & b_{n-1} & c_{n-1} & d_{n-1} \
   & & & & & a_{n} & b_{n} & c_{n}
   \end{bmatrix}
   $$

本文将使用matlab(R2022a)（win11 64位操作系统）逐一实现上列线性方程组的追赶法。详细原理及算法分析来源于最后的参考资料，文中并不提及。

三对角

算法步骤
一、替换(角标表示替换后的元素，下文同理)
$$
c^{}{i} = \begin{cases}
\frac{c{i}}{b_{i}} & \text{if }i = 1, \
\frac{c_{i}}{b_{i}-a_{i}c^{}_{i}} & \text{if }i = 2,3,\ldots,n-1.
\end{cases}
$$
$$
d^{}{i} =\begin{cases}
\frac{d{i}}{b_{i}} & \text{if }i = 1, \
\frac{d_{i}-a_{i}d^{}{i-1}}{b{i}-a_{i}c^{}{i-1}} & \text{if }i = 2,3,\ldots,n.
\end{cases}
$$
其中$$d{i} \in y$$.

二、计算解向量（追与赶）
$$
x_{n} = d^{}{n}, x{i} = d^{}{i}-c{i}x_{i+1} \text{ for } i = n-1,n-2,\ldots,1.
$$

代码实现

Tridiagonal.m

function x = Tridiagonal(a, b, c, d)
% a、b、c 分别为三对角线的下对角线、主对角线和上对角线
% d 为常数向量
n = length (d); % 矩阵维度
% 追赶法的第一步：消元
c(1) = c(1)/b(1);
d(1) = d(1)/b(1);
for i = 2 : n-1
    c(i) = c(i) / (b(i)-a(i)*c(i-1));
    d(i) = (d(i)-a(i-1)*d(i-1)) / (b(i)-a(i-1)*c(i-1));
end
d(n) = (d(n)-a(n-1)*d(n-1)) / (b(n)-a(n-1)*c(n-1));
% 追赶法的第二步：回代
x(n) = d (n);
for i = n-1:-1:1
    x(i) = d(i)-c(i)*x(i+1);
end
end

循环Toeplitz三对角

算法步骤
一、$$A = LU$$, 其中

$$
L = \begin{bmatrix}
d_1 & & & & & & \
a_2 & d_2 & & & & & \
& \ddots & \ddots & & & & \
& & a_{n-2} & d_{n-2} & & & \
& & & a_{n-1} & d_{n-1} & & \
\alpha_1 & \alpha_2 & \cdots & \alpha_{n-2} & \alpha_{n-1} & d_n &
\end{bmatrix}, \
U = \begin{bmatrix}
1 & u_1 & & & & \beta_1 \
& 1 & u_2 & & & \beta_2 \
& & \ddots & \ddots & & \vdots \
& & & 1 & u_{n-2} & \beta_{n-2} \
& & & & 1 &\beta_{n-1}\
& & & & & 1
\end{bmatrix}
$$
其中
$$
d_1 = b_1, \quad u_1 = c_1/d_1, \quad \alpha_1 = c_n, \quad \beta_1 = \alpha_1/d_1,
$$
$$
\begin{cases}
d_i &= b_i - a_i u_{i-1} \
u_i &= c_i/d_i , \quad i = 2,3,\dots,n-2 \
\alpha_i &= -\alpha_{i-1}u_{i-1}
\end{cases} \
$$
$$
d_{n-1} = b_{n-1} - \alpha_{n-2}u_{n-2}, \
\alpha_{n-1} = a_n - \alpha_{n-2}u_{n-2}, \
\beta_{n-1} = (c_{n-1} - a_{n-1}\beta_{n-2})/d_{n-1}, \
d_n = b_n - \sum_{i=1}^{n-1}\alpha_i\beta_i.
$$
二、首先求解方程$$Ly=f$$
$$
\begin{aligned}
\begin{cases}
y_1 &= f_1/d_1, \
y_i &= (f_i - a_i y_{i-1})/d_i, \quad i = 2,3,\dots,n-1, \
y_n &= (f_n - \sum_{i=1}^{n-1}\alpha_i y_i)/d_n.
\end{cases}
\end{aligned}
$$
三、然后求解方程$$Ux=y$$
$$
\begin{aligned}
\begin{cases}
x_n &= y_n, \
x_{n-1} &= y_{n-1} - \beta_{n-1}x_n, \
x_i &= y_i - u_ix_{i+1} - \beta_ix_n, \quad i = n-2,n-3,\dots,2,1.
\end{cases}
\end{aligned}
$$

代码实现

CTTChase.m

function x = CTTChase(a,b,c,f)
n = length (f); % 矩阵维度
c(1) = c(1)/b(1);
alpha = zeros(n-1,1);
alpha(1) = c(n); 
beta = zeros(n-1,1);
beta(1) = alpha(1)/b(1);
for i = 2:n-2
    b(i) = b(i) - a(i)*c(i-1);
    c(i) = c(i)/b(i);
    alpha(i) = - alpha(i-1)*c(i-1);
    beta(i) = - (beta(i-1)*a(i))/b(i);
end
b(n-1) = b(n-1) - alpha(n-2)*c(n-2);
alpha(n-1) = a(n) - alpha(n-2)*c(n-2);
beta(n-1) = (c(n-1)-a(n-1)*beta(n-2))/b(n-1);
b(n) = b(n) - sum(alpha.*beta);

f(1) = f(1)/b(1);
for i = 2:n-1
    f(i) = (f(i)-a(i)*f(i-1))/b(i);
end
sum1 = 0;
for i = 1:n-1
    sum1 = sum1 + alpha(i)*f(i);
end
f(n) = (f(n)-sum1)/b(n);

x(n) = f(n);
x(n-1) = f(n) - beta(n-1)*x(n);
for i = n-2:-1:1
    x(i) = f(i) - c(i)*x(i+1) - beta(i)*x(n);
end
end

五对角

算法步骤
一、替换
$$
b^{}{i} =
\begin{cases}
b{i} & \text{if }i = 2, \
b_{i} - a_{i}d^{}{i-2} & \text{if }i = 3,4,\ldots,n.
\end{cases} \
$$
$$
c^{*}{i} =
\begin{cases}
c_{i} & \text{if }i = 1, \
c_{i} - \frac{b^{}{i}d{i-1}}{c_{i-1}} & \text{if }i = 2, \
c_{i} - a_{i}e^{}{i-2} - b^{*}{i}d^{}_{i-1} & \text{if }i = 3,4,\ldots,n.
\end{cases} \
$$
$$
d^{}{i} =
\begin{cases}
\frac{d{i}}{c_{i}} & \text{if }i = 1, \
\frac{d_{i} - b^{}_{i}d^{}{i-1}}{c^{*}{i}} & \text{if }i = 2,3,\ldots,n-1.
\end{cases} \
$$
$$
e^{}{i} =
\frac{e{i}}{c^{}_{i}} \quad i = 1,2,\ldots,n-2.
$$

二、设 $$ Ly = f $$，求$$ y $$.
$$
y^{}{i} =
\begin{cases}
\frac{y{i}}{c_{i}} & \text{if }i = 1, \
\frac{y_{i} - b^{}{i}y{i-1}}{c^{}{i}} & \text{if }i = 2 \
\frac{y{i} - a_{i}y^{}{i-2} - b^{*}{i}y^{}_{i-1}}{c^{}{i}} & \text{if }i = 3,4,\ldots,n.
\end{cases}
$$
三、计算解向量$x$
$$
\begin{cases}
x{n} = y^{}{n} \
x{n-1} = y^{}{n-1} - d^{*}{n-1}x_{n} \
x_{i} = y^{}_{i} - d^{}{i}x{i+1} - e^{*}{i}x{i+2} \text{ for } i = n-2,n-3,\ldots,1.
\end{cases}
$$
代码实现

Pentadiagonal.m

function x = Pentadiagonal(a, b, c, d, e, f)
a=[0,0,a];
b=[0,b];
d=[d,0];
e=[e,0,0];
n = length (f); % 矩阵维度
alpha=zeros(n,1); %1~n
ganma=zeros(n,1);%2~n
z=zeros(n,1);%3~n

beta=zeros(n,1);%1~n-1
q=zeros(n,1);%1~n-2

alpha(1)=c(1);
beta(1)=d(1)/alpha(1);
% beta(1)=d(1)/d(2);

ganma(2)=b(2);
alpha(2)=c(2)-ganma(2)*beta(1);

for i=1:2
    q(i)=e(i)/alpha(i);
end

beta(2)=(d(2)-ganma(2)*q(1))/alpha(2);

for i=3:n
    if i<=n-2
        z(i)=a(i);
        ganma(i)=b(i)-z(i)*beta(i-2);
        alpha(i)=c(i)-z(i)*q(i-2)-ganma(i)*beta(i-1);
        q(i)=e(i)/alpha(i);
        beta(i)=(d(i)-ganma(i)*q(i-1))/alpha(i);
    elseif i<=n-1
        z(i)=a(i);
        ganma(i)=b(i)-z(i)*beta(i-2);
        alpha(i)=c(i)-z(i)*q(i-2)-ganma(i)*beta(i-1);
        beta(i)=(d(i)-ganma(i)*q(i-1))/alpha(i);
    else
        z(i)=a(i);
        ganma(i)=b(i)-z(i)*beta(i-2);
        alpha(i)=c(i)-z(i)*q(i-2)-ganma(i)*beta(i-1);
    end
end

y=zeros(n,1);
y(1)=f(1)/alpha(1);
y(2)=(f(2)-ganma(2)*y(1))/alpha(2);
for i=3:n
    y(i)=(f(i)-z(i)*y(i-2)-ganma(i)*y(i-1))/alpha(i);
end

x=zeros(n,1);
x(n)=y(n);
x(n-1)=y(n-1)-beta(n-1)*x(n);
for i=n-2:-1:1
    x(i)=y(i)-q(i)*x(i+2)-beta(i)*x(i+1);
end
end

Latex公式源码文件

下载地址：点我

参考资料

[1] 金一庆、陈越、王冬梅，数值方法，机械工业出版社，2000
[2] 李青、王能超，解循环三对角线性方程组的追赶法，小型微型计算机系统，第23卷第11期，2002年11月
[3] 王礼广，蔡放，熊岳山，五对角线性方程组追赶法，南华大学学报(自然科学版)，第22卷第1期，2008年3月