(追赶法求三对角矩阵、LU分解)

(追赶法求三对角矩阵、LU分解)

应用场景

数模中的微分方程问题中偶尔回涉及到使用向后差分法求解问题，不论是六点法向后差分还是四点法向后差分都会涉及到求解三对角矩阵的情况。

以四点法向后差分为例

四点法得到的代数矩阵为三对角线性方程组AU=b，一般可以用matlab求逆命令求解U，但因本题中系数矩阵规模较大，且主对角线和两个次对角线之外，其余元素均为0，因而求逆解方程组会导致浪费内存、程序运行缓慢。因此本文利用追赶法求解三对角线性方程组。三对角线性方程组的追赶法求解过程如下图。

计算过程

三对角矩阵的形式如下：

需要将三对角阵分解成如下的形式：

求解Ax=d等价于求解Ly=d和Ux=y.

因而经过推导得到解三对角线性方程组的追赶法公式.

matlab代码如下

% x为待求系数向量，d为右端已知常量向量，长度为n
% a为-1对角向量，b为对角向量，c为1对角向量，ac长度为n-1，b长度为n 
function x=machase(d,n,a,b,c)
y=zeros(n,1);
x=zeros(n,1);
alpha=ones(n,1);
beta=ones(n-1,1);
alpha(1)=b(1);beta(1)=c(1)/b(1);
for i=2:n-1
    alpha(i)=b(i)-a(i-1)*beta(i-1);
    beta(i)=c(i)/alpha(i);
end
alpha(n)=b(n)-a(n-1)*beta(n-1);
y(1)=d(1)/b(1);
for i=2:n
    y(i)=(d(i)-a(i-1)*y(i-1))/alpha(i);
end
x(n)=y(n);
for i=(n-1):-1:1 
    x(i)=y(i)-beta(i)*x(i+1);
end
end

测试命令

a=-1*ones(4,1);
b=2*ones(5,1);
d=[1 0 0 0 0];
machase(d,5,a,b,a)

结果

    0.8333
    0.6667
    0.5000
    0.3333
    0.1667