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

追赶法求解三对角线性方程组
一 实验目的
利用编程方法实现追赶法求解三对角线性方程组。

二 实验内容
1、 学习和理解追赶法求解三对角线性方程组的原理及方法；
2、 利用MATLAB 编程实现追赶法；
3、 举例进行求解，并对结果进行分。

三 实验原理
设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵
11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦
………… 这种方程组称为三对角线性方程组。

显然，A 是上下半宽带都是1的带状矩阵。

设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。

其中L 是下三角矩阵，U 是单位上三角矩阵。

利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………
得到：
由上列各式可以得到L 和U 。

引入中间量y ，令
y Ux =，则有：
已知
L 和d ，可求得y 。

则可得到y 的求解表达式：
11/1
2,3,,()(1)*y()=()[()(1)]/y d l i n
m i y i li i di y i di m i y i li
==-+=--…
1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui
=*===≤≤+=+++≤≤-=•=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
……………
由y Ux =得：
111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
………… 可得到X 的求解表达式：
()()
1,2,,1
()()u()(1)
x n y n i n n x i y i i x i ==--=-+…
从而得到Ax=d 的解x 。

四 Matlab 编程
根据以上的实验原理，在Matlab 中编程如下函数x=Trid （A,d ）
在Matlab中新建Trid.m,其中程序为如上虚线框内代码，放在工
作目录。

在Command Window输入以下语句：
clear all;clc;
fprintf('输入非奇异三对角系数矩阵A\n');
A=input('A matrix='); % 输入系数矩阵
fprintf('系数矩阵');A
if det(A)==0 % 判断系数矩阵是否奇异
fprintf('系数矩阵A奇异！！！请重新输入！\n');
fprintf('重新输入非奇异三对角系数矩阵A\n');
A=input('A matrix=');
fprintf('系数矩阵');A
end
fprintf('输入矩阵d\n');
d=input('d matrix=');
fprintf('矩阵d');d
x=Trid(A,d) % 调用Trid.m中的Trid函数进行求解
fid = fopen('Ax=d.txt', 'wt'); % 生成Ax=d.txt文件
fprintf(fid,'%s\r\n','利用三对角线追赶法求解Ax=d');
fprintf(fid,'%s\r\n','==================================');
for i=1:size(A,1)
fprintf(fid, '%.1f\t', A(i,:)); % 输出Ax=d，以上=为分隔符fprintf(fid, '%s', 'x');
fprintf(fid, '%u\t', i);
fprintf(fid, '%.1f\n', d(i));
end
fprintf(fid,'%s\r\n','==================================');
五 举例计算及分析
以课本（《数值分析》第4版，颜庆津，北京航空航天大学出版社）27页例3为例进行计算，输入系数矩阵A 和d ：
41141=4114A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
1411，10.5=132d ⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦ 调用x=Trid(A,d)后，并生成Ax=d.txt 文件，Command Window
同时也会输出解为0.20.2=0.50.80.3x ⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦
，与课本答案一致; Ax=d.txt 文件内容
如下图：
fprintf(fid,'%s\r\n','求解得到结果如下：'); for i=1:size(A,1)
fprintf(fid,'%s','x'); % 输出解向量x （i ）
fprintf(fid, '%u', i);
fprintf(fid, '%s\t', '='); fprintf(fid,'%.5f\r\n',x(i)); end fclose(fid)
再次对程序进行验证，输入矩阵如下：
12213=120.531A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦，21=13d ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦
利用matlab 计算Ax 得到Ax=d ，验证所得到的x 即为方程组的解。

由以上两组计算可表明，该程序能满足追赶法解三对角线性方程
组，其中要求系数矩阵满足要求，即非奇异，三对角且前（n-1）个顺序主子式都不为零。

在又换了一组系数矩阵后，系数矩阵如下：
A=[1 1 0 0 0;2 2 1 0 0;0 2 3 1 0;0 0 2 4 1;0 0 0 2 5]
即11221=23124125A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
，11=111d ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ 运行程序后出错，查找问题发现A 进行Crout 分解时，
l1=1
m2=2
u1=1
l2=a2-m2*u1=2-2×1=0
再计算u2=c2/l2时无法计算，分析发现A 的2阶顺序主子式为零，无法进行Crout 分解，所以该程序还无法判断输入的系数矩阵的（n-1）阶顺序主子式是否为零，这是不足的地方，也无法判断系数矩阵是否是三对角矩阵，还需要在程序上进行修改完善。

六 实验总结
1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………
通过本次的课程实验设计，对追赶法求解三对角线性方程组有了深刻的认识和了解，用matlab编程实现了该解法；同时在实验过程中明白了系数矩阵需要满足一定的要求才能使用该种追赶法方法。