解线性方程组——追赶法解三对角方程组 | 北太天元

一、问题描述

对于线性方程组
$$Ax=b,A=\left(\begin{array}{cccccc}{b}_1& c_1& & & & \\ a_2& b_2& c_2& & & \\ & \ddots & \ddots & \ddots & & \\ & & \ddots & \ddots & \ddots & \\ & & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & & a_n& b_n\end{array}\right),b=\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right)$$

其中
$$\{\begin{array}{l}\begin{array}{rl}& |b_1|>|c_1|>0\\ & |b_n|>|a_n|>0\\ & |b_i|\ge |a_i|+|c_i|,a_i{c}_i\ne 0,i=2,3,\cdots n-1\end{array}\end{array}$$

二、追赶法解三对角方程组

$$\begin{array}{c}A=\left(\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right),L=\left(\begin{array}{ccccc}1& & & & \\ l_2& 1& & & \\ & l_3& \ddots & & \\ & & \ddots & 1& \\ & & & l_n& 1\end{array}\right),U=\left(\begin{array}{ccccc}{u}_1& d_1& & & \\ & u_2& d_2& & \\ & & \ddots & \ddots & \\ & & & u_{n-1}& d_{n-1}\\ & & & & u_n\end{array}\right)\end{array}$$

$$A=LU$$

首先通过前两行与前两列可以知道
$$b_1=u_1,a_2=l_2{u}_1,c_1=d_1$$得 $u_1,d_1,l_2$

再从下式看两边是如何对应上的

$$A:\begin{array}{cccc}& i-1& i& i+1\\ i& a_i& b_i& c_i\end{array}$$

$$L:\begin{array}{cccc}& i-1& i& i+1\\ i& l_i& 1& \end{array}U:\begin{array}{cccc}& i-1& i& i+1\\ i-2& d_{i-2}& & \\ i-1& u_{i-1}& d_{i-1}& \\ i& & u_i& d_i\\ i+1& & & u_{i+1}\end{array}$$

可以看出
$$\begin{array}{rl}& a_i=l_i{u}_{i-1}\\ & b_i=l_i{d}_{i-1}+u_i\\ & c_i=d_i\end{array}$$

$d_i$已经得出，接下来只需算 $l_i$ 与 $u_i$ 即可

再有$u_1=b_1$, $i=2\to n$

$$\begin{array}{c}{l}_i=\frac{a_i}{u_{i-1}}\\ u_i=b_i-l_i\cdot c_{i-1}\end{array}$$

这样子, $L$ 和 $U$就得到了

接下来，解 $Ly=b$ , $Ux=y$,这一步可以用之前写好的程序直接解出

由于这里，结果比较简易，直接将结果写了出来：

解 $Ly=b$： （追）
$$\left(\begin{array}{ccccc}1& & & & \\ l_2& 1& & & \\ & l_3& \ddots & & \\ & & \ddots & 1& \\ & & & l_n& 1\end{array}\right)(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_{n-1}\\ y_n\end{array})=\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_{n-1}\\ f_n\end{array}\right)$$
$$y_1=f_1,y_i=f_i-l_i\cdot y_{i-1}(i=2,3\cdots n)$$

解 $Ux=y$：（赶）
$$\left(\begin{array}{ccccc}{u}_1& d_1& & & \\ & u_2& d_2& & \\ & & \ddots & \ddots & \\ & & & u_{n-1}& d_{n-1}\\ & & & & u_n\end{array}\right)(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n-1}\\ x_n\end{array})=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_{n-1}\\ y_n\end{array}\right)$$
$$x_n=\frac{y_n}{u_n},x_i=\frac{(y_i-d_i\cdot x_{i+1})}{u_i}(i=n-1,\cdots ,2,1)$$

这样 $x$ 就求出来了

三、算法

四、北太天元源代码

function [x]=tridiag_chase(A,f)
    % 追赶法解三对角方程组
    % 输入:  适用的三对角矩阵A,  右端向量f
    % 输出:  解，列向量的形式 x
    % 创建时间: 1/26/2024
    % 版本: 1.0
    n = length(A);
    % 将三对角提取出来
    b = diag(A,0); a = diag(A,-1); c = diag(A,1);
    % 处理一下角标 a是从a_2开始, l从l_2 开始
    a = cat(1,[0],a); % a = [0, diag(A,-1)]
    u = zeros(1,n);l = zeros(1,n);
    u(1) = b(1);
    for i = 2:1:n
        l(i) = a(i)/u(i-1);
        u(i) = b(i)-l(i)*c(i-1);
    end
    % Ly = b
    y = zeros(1,n);
    y(1) = f(1);
    for i =2:1:n
        y(i) = f(i)-l(i)*y(i-1);
    end
    % Ux= y
    x = zeros(n,1);
    x(n) = y(n)/u(n);
    for i =n-1:-1:1
        x(i) = (y(i) - c(i)* x(i+1))/u(i);
    end
end

保存为 tridiag_chase.m 文件

简单测试一下

在 10、100、500 、1000阶 三对角方程组下, 追赶法与gauss列主元消去法 的时间消耗对比

% file need: tridiag_chase.m,  gsem_column.m
% time: 1/26/2024
%%
clc,clear all;
n = 10;  % 方程组的阶数 10 100 500 1000
A = diag(2*ones(n,1)) + diag(-1*ones(n-1,1),1) + diag(-1*ones(n-1,1),-1);
f = (1:n)';
t1 = tic;
x1 = tridiag_chase(A,f);
toc(t1);

t2 = tic;
x2 = gsem_column(A,f);
toc(t2);

delta = norm(abs(x1-x2));

n	追赶法	Gauss列主元	delta
10	0.242069 s	0.253921 s	0
100	0.249970 s	0.997883 s	0
500	0.373002 s	19.068616 s	0
1000	0.440046 s	81.370982 s	0

随着方程组阶数的增大、其优势愈加明显。

其中用到的 Gauss列主元消去法 ，见左方蓝色字体。