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

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

一 实验目的

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

二 实验内容

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)