【matlab】【数值分析】针对特殊矩阵的追赶法的matlab实现

原文链接：
点我！这就是本人的博客喵，快来看喵！

下面的追赶法算法原理不予介绍，在参考文献中有原文可参照。

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

1. 常对角矩阵(Toeplitz矩阵)
   $$A=\left[\begin{array}{cccccc}{a}_0& a_{-1}& a_{-2}& \dots & \dots & a_{-n+1}\\ a_1& a_0& a_{-1}& \ddots & & ⋮\\ a_2& a_1& \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & a_{-1}& a_{-2}\\ ⋮& & \ddots & a_1& a_0& a_{-1}\\ a_{n-1}& \dots & \dots & a_2& a_1& a_0\end{array}\right]$$
2. 循环矩阵
   $$A=\left[\begin{array}{cccccc}{a}_0& a_{n-1}& a_{n-2}& \dots & \dots & a_1\\ a_1& a_0& a_{n-1}& a_{n-2}& & ⋮\\ a_2& a_1& a_0& a_{n-1}& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & a_1& a_{n-2}\\ ⋮& & \ddots & a_1& a_0& a_{n-1}\\ a_{n-1}& a_{n-2}& a_{n-3}& \dots & \dots & a_0\end{array}\right]$$
   为特殊的Toeplitz矩阵。
3. 三对角矩阵
   $$A=\left[\begin{array}{ccccc}{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{array}\right]$$
4. 循环Toeplitz三对角线性方程组
   $$A=\left[\begin{array}{ccccc}{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{array}\right]$$
5. 五对角矩阵
   $$A=\left[\begin{array}{cccccccc}{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{array}\right]$$

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

三对角

算法步骤
一、替换(*角标表示替换后的元素，下文同理)
$$c_i^{\ast }=\{\begin{array}{ll}\frac{c_i}{b_i}& \text{if }i=1,\\ \frac{c_i}{b_i-a_i{c}_i^{\ast }}& \text{if }i=2,3,\dots ,n-1.\end{array}$$
$$d_i^{\ast }=\{\begin{array}{ll}\frac{d_i}{b_i}& \text{if }i=1,\\ \frac{d_i-a_i{d}_{i-1}^{\ast }}{b_i-a_i{c}_{i-1}^{\ast }}& \text{if }i=2,3,\dots ,n.\end{array}$$
其中$$d_i\in y$$.

二、计算解向量（追与赶）
$$x_n=d_n^{\ast },x_i=d_i^{\ast }-c_i{x}_{i+1}\text{ for }i=n-1,n-2,\dots ,1.$$

代码实现

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$$, 其中

$$\begin{gathered} L=\left[\begin{array}{ccccccc}{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{array}\right], \\ U=\left[\begin{array}{cccccc}1& u_1& & & & {\beta }_1\\ & 1& u_2& & & {\beta }_2\\ & & \ddots & \ddots & & ⋮\\ & & & 1& u_{n-2}& {\beta }_{n-2}\\ & & & & 1& {\beta }_{n-1}\\ & & & & & 1\end{array}\right] \end{gathered}$$
其中
$$d_1=b_1,u_1=c_1/d_1,{\alpha }_1=c_n,{\beta }_1={\alpha }_1/d_1,$$
$$\{\begin{array}{ll}{d}_i& =b_i-a_i{u}_{i-1}\\ u_i& =c_i/d_i,i=2,3,\dots ,n-2\\ {\alpha }_i& =-{\alpha }_{i-1}{u}_{i-1}\end{array}$$
$$\begin{gathered} 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. \end{gathered}$$
二、首先求解方程$$Ly=f$$
$$\begin{array}{r}\{\begin{array}{ll}{y}_1& =f_1/d_1,\\ y_i& =(f_i-a_i{y}_{i-1})/d_i,i=2,3,\dots ,n-1,\\ y_n& =(f_n-{\sum }_{i=1}^{n-1}{\alpha }_i{y}_i)/d_n.\end{array}\end{array}$$
三、然后求解方程$$Ux=y$$
$$\begin{array}{r}\{\begin{array}{ll}{x}_n& =y_n,\\ x_{n-1}& =y_{n-1}-{\beta }_{n-1}{x}_n,\\ x_i& =y_i-u_i{x}_{i+1}-{\beta }_i{x}_n,i=n-2,n-3,\dots ,2,1.\end{array}\end{array}$$

代码实现

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^{\ast }=\{\begin{array}{ll}{b}_i& \text{if }i=2,\\ b_i-a_i{d}_{i-2}^{\ast }& \text{if }i=3,4,\dots ,n.\end{array}$$
$$c_i^{\ast }=\{\begin{array}{ll}{c}_i& \text{if }i=1,\\ c_i-\frac{b_i^{\ast }{d}_{i-1}}{c_{i-1}}& \text{if }i=2,\\ c_i-a_i{e}_{i-2}^{\ast }-b_i^{\ast }{d}_{i-1}^{\ast }& \text{if }i=3,4,\dots ,n.\end{array}$$
$$d_i^{\ast }=\{\begin{array}{ll}\frac{d_i}{c_i}& \text{if }i=1,\\ \frac{d_i-b_i^{\ast }{d}_{i-1}^{\ast }}{c_i^{\ast }}& \text{if }i=2,3,\dots ,n-1.\end{array}$$
$$e_i^{\ast }=\frac{e_i}{c_i^{\ast }}i=1,2,\dots ,n-2.$$

二、设 $$Ly=f$$，求$$y$$.
$$y_i^{\ast }=\{\begin{array}{ll}\frac{y_i}{c_i}& \text{if }i=1,\\ \frac{y_i-b_i^{\ast }{y}_{i-1}}{c_i^{\ast }}& \text{if }i=2\\ \frac{y_i-a_i{y}_{i-2}^{\ast }-b_i^{\ast }{y}_{i-1}^{\ast }}{c_i^{\ast }}& \text{if }i=3,4,\dots ,n.\end{array}$$
三、计算解向量$x$
$$\{\begin{array}{l}{x}_n=y_n^{\ast }\\ x_{n-1}=y_{n-1}^{\ast }-d_{n-1}^{\ast }{x}_n\\ x_i=y_i^{\ast }-d_i^{\ast }{x}_{i+1}-e_i^{\ast }{x}_{i+2}\text{ for }i=n-2,n-3,\dots ,1.\end{array}$$
代码实现

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月