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

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

一 实验目的

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

二 实验内容

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 ）

n=size(A,1); % n为系数矩阵的行数

a(1)=A(1,1)

b(1)=0

c(1)=A(1,2)

for i=2:n-1

a(i)=A(i,i)

b(i)=A(i,i-1);

c(i)=A(i,i+1);

end

a(n)=A(n,n)

b(n)=A(n,n-1)

c(n)=0

l(1)=a(1); %开始求解L，U

m(1)=0

for i=2:n

m(i)=b(i) %求得m(i)

u(i-1)=c(i-1)/l(i-1); %求得u(i)

l(i)=a(i)-b(i)*u(i-1); %求得l(i)

end

u(n)=0

y(1)=d(1)/l(1);

for i=2:n

y(i)=[d(i)-m(i)*y(i-1)]/l(i); %求得y(i) end

x(n)=y(n);

for i=n-1:-1:1

x(i)=y(i)-u(i)*x(i+1); %求得x(i) end

x=x' %将x转置，变为列向量

在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','==================================');