解线性方程组——追赶法解三对角方程组

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

一、问题描述

对于线性方程组
$Ax=b,\quad A=\begin{pmatrix}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{pmatrix},\quad b=\begin{pmatrix} f_1\\f_2\\ \vdots\\ f_n \end{pmatrix}$

其中
$\begin{cases} \begin{aligned}&\quad|b_1|>|c_1|>0\\&\quad|b_n|>|a_n|>0\\&\quad|b_i|\geq|a_i|+|c_i|,\quad a_ic_i\neq0,i=2,3,\cdots n-1\end{aligned} \end{cases}$

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

$\begin{equation*} A=\begin{pmatrix}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{pmatrix}, {L}=\begin{pmatrix}1\\l_2&1\\&l_3&\ddots\\&&\ddots&1\\&&&l_n&1\end{pmatrix},{U}=\begin{pmatrix}u_1&d_1\\&u_2&d_2\\&&\ddots&\ddots\\&&&u_{n-1}&d_{n-1}\\&&&&u_n\end{pmatrix} \end{equation*}$

$A = LU$

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

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

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

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

可以看出
$\begin{aligned} &a_i = l_i u_{i-1}\\ &b_i = l_i d_{i-1} +u_i\\ &c_i = d_i \end{aligned}$

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

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

$\begin{gathered} l_i = \dfrac{a_i}{u_{i-1}}\\ u_i = b_i - l_i \cdot c_{i-1} \end{gathered}$

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

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

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

解 $L y = b$ ：（追）
$\left.\begin{pmatrix}1\\l_2&1\\&l_3&\ddots\\&&\ddots&1\\&&&l_n&1\end{pmatrix}\left(\begin{array}{c}y_1\\y_2\\\vdots\\y_{n-1}\\y_n\end{array}\right.\right)=\left(\begin{array}{c}f_1\\f_2\\\vdots\\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{pmatrix}u_1&d_1&&&\\&u_2&d_2&&\\&&\ddots&\ddots&\\&&&u_{n-1}&d_{n-1}\\&&&&u_n\end{pmatrix}\left(\begin{array}{c}x_1\\x_2\\\vdots\\x_{n-1}\\x_n\end{array}\right.\right)=\left(\begin{array}{c}y_1\\y_2\\\vdots\\y_{n-1}\\y_n\end{array}\right)$
$x_n = \dfrac{y_n}{u_n}, x_i = \dfrac{(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.0n = 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:nl(i) = a(i)/u(i-1);u(i) = b(i)-l(i)*c(i-1);end% Ly = by = zeros(1,n);y(1) = f(1);for i =2:1:ny(i) = f(i)-l(i)*y(i-1);end% Ux= yx = zeros(n,1);x(n) = y(n)/u(n);for i =n-1:-1:1x(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));