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

追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法；2、 利用MATLAB编程实现追赶法；3、 举例进行求解，并对结果进行分。三 实验原理设n元线性方程组的系数矩阵A为非奇异的三对角矩阵这种方程组称为三对角线性方程组。显然，A是上下半宽带都是1的带状矩阵。设A的前n-1个顺序主子式都不为零，根据定理2.5的推论，A有唯一的Crout分解，并且是保留带宽的。其中L是下三角矩阵，U是单位上三角矩阵。利用矩阵相乘法，可以得到：由上列各式可以得到L和U。 引入中间量y，令，则有： 已知L和d，可求得y。则可得到y的求解表达式：由得：可得到X的求解表达式：从而得到的解x。四 Matlab编程 function x=Trid(A,d)% 追赶法求解三对角的线性方程组 Ax=d% b为主对角线元素，a，c分别为次对角线元素，d为右端项 % A= a1 c1% b2 a2 c2% .% b(n-1) a(n-1) c(n-1)% b(n) a(n) % a=a1.a(n) % 把系数矩阵的三对角转变成3个列向量 % b=0 b2.b(n) % 不足的元素用0代替% c=c1.c(n-1) 0根据以上的实验原理，在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);enda(n)=A(n,n)b(n)=A(n,n-1)c(n)=0l(1)=a(1); %开始求解L，Um(1)=0for 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)endu(n)=0y(1)=d(1)/l(1);for i=2:n y(i)=d(i)-m(i)*y(i-1)/l(i); %求得y(i)endx(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(输入非奇异三对角系数矩阵An);A=input(A matrix=); % 输入系数矩阵 fprintf(系数矩阵);Aif det(A)=0 % 判断系数矩阵是否奇异fprintf(系数矩阵A奇异！请重新输入！n);fprintf(重新输入非奇异三对角系数矩阵An);A=input(A matrix=); fprintf(系数矩阵);Aend fprintf(输入矩阵dn);d=input(d matrix=);fprintf(矩阵d);dx=Trid(A,d) % 调用Trid.m中的Trid函数进行求解fid = fopen(Ax=d.txt, wt); % 生成Ax=d.txt文件fprintf(fid,%srn,利用三对角线追赶法求解Ax=d);fprintf(fid,%srn,=);for i=1:size(A,1)fprintf(fid, %.1ft, A(i,:); % 输出Ax=d，以上=为分隔符fprintf(fid, %s, x);fprintf(fid, %ut, i);fprintf(fid, %.1fn, d(i);end fprintf(fid,%srn,=);fprintf(fid,%srn,求解得到结果如下：);for i=1:size(A,1)fprintf(fid,%s,x); % 输出解向量x（i）fprintf(fid, %u, i);fprintf(fid, %st, =);fprintf(fid,%.5frn,x(i);endfclose(fid)五 举例计算及分析以课本（数值分析第4版，颜庆津，北京航空航天大学出版社）27页例3为例进行计算，输入系数矩阵A和d：，调用x=Trid(A,d)后，并生成Ax=d.txt文件，Command Window同时也会输出解为，与课本答案一致; Ax=d.txt文件内容如下图： 再次对程序进行验证，输入矩阵如下：， 利用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即，运行程序后出错，查找问题发现A进行Crout分解时，l1=1m2=2u1=1l2=a2-m2*u1=2-21=0再计算u2=c2/l2时无法计算，分析发现A的2阶顺序主子式为零，无法进行Crout分解，所以该程序还无法判断输入的系数矩阵的（n-1）阶顺序主子式是否为零，这是不足的地方，也无法判断系数矩阵是否