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

《数值分析》数值实验XXX0追赶法求解三对角线性方程组一实验目的利用编程方法实现追赶法求解三对角线性方程组。二实验内容1、学习和理解追赶法求解三对角线性方程组的原理及方法；2、利用MATLAB编程实现追赶法；3、举例进行求解，并对结果进行分。三实验原理设n元线性方程组Ax=d的系数矩阵A为非奇异的三对角矩阵11222=(1)(n1)()()acbacAancbnan…………这种方程组称为三对角线性方程组。显然，A是上下半宽带都是1的带状矩阵。设A的前n-1个顺序主子式都不为零，根据定理2.5的推论，A有唯一的Crout分解，并且是保留带宽的。其中L是下三角矩阵，U是单位上三角矩阵。利用矩阵相乘法，可以1112212(1)1u(n1)()()1lumluALUlnmnln……………《数值分析》数值实验XXX1得到：由上列各式可以得到L和U。引入中间量y，令yUx，则有：已知L和d，可求得y。则可得到y的求解表达式：11/12,3,,()(1)*y()=()[()(1)]/ydlinmiyiliidiyidimiyili…1111111/1(2)(1)(1)u(1)(11)/(1)(1)(1)lalucuclmibiinaimiiliinciliuiuicililiaibiuiAxLUxLydLyd1112222(1)(n1)(n1)()()(n)(n)lydmlydlnydmnlnyd……………《数值分析》数值实验XXX2由yUx得：111112221u(n1)(n1)(n1)1(n)(n)uxyuxyxyxy…………可得到X的求解表达式：()()1,2,,1()()u()(1)xnyninnxiyiixi…从而得到Ax=d的解x。四Matlab编程根据以上的实验原理，在Matlab中编程如下函数x=Trid（A,d）functionx=Trid(A,d)%追赶法求解三对角的线性方程组Ax=d%b为主对角线元素，a，c分别为次对角线元素，d为右端项%A=[a1c1%b2a2c2%......%b(n-1)a(n-1)c(n-1)%b(n)a(n)]%a=[a1...a(n)]%把系数矩阵的三对角转变成3个列向量%b=[0b2...b(n)]%不足的元素用0代替%c=[c1...c(n-1)0]《数值分析》数值实验XXX3n=size(A,1);%n为系数矩阵的行数a(1)=A(1,1)b(1)=0c(1)=A(1,2)fori=2:n-1a(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)=0fori=2:nm(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);fori=2:ny(i)=[d(i)-m(i)*y(i-1)]/l(i);%求得y(i)endx(n)=y(n);fori=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);%求得x(i)endx=x'%将x转置，变为列向量《数值分析》数值实验XXX4在Matlab中新建Trid.m,其中程序为如上虚线框内代码，放在工作目录。在CommandWindow输入以下语句：clearall;clc;fprintf('输入非奇异三对角系数矩阵A\n');A=input('Amatrix=');%输入系数矩阵fprintf('系数矩阵');Aifdet(A)==0%判断系数矩阵是否奇异fprintf('系数矩阵A奇异！！！请重新输入！\n');fprintf('重新输入非奇异三对角系数矩阵A\n');A=input('Amatrix=');fprintf('系数矩阵');Aendfprintf('输入矩阵d\n');d=input('dmatrix=');fprintf('矩阵d');dx=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','==================================');fori=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));endfprintf(fid,'%s\r\n','==================================');《数值分析》数值实验XXX5五举例计算及分析以课本（《数值分析》第4版，颜庆津，北京航空航天大学出版社）27页例3为例进行计算，输入系数矩阵A和d：41141=4114A1411，10.5=132d调用x=Trid(A,d)后，并生成Ax=d.txt文件，CommandWindow同时也会输出解为0.20.2=0.50.80.3x，与课本答案一致;Ax=d.txt文件内容如下图：fprintf(fid,'%s\r\n','求解得到结果如下：');fori=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));endfclose(fid)《数值分析》数值实验XXX6再次对程序进行验证，输入矩阵如下：12213=120.531A，21=13d利用matlab计算Ax得到Ax=d，验证所得到的x即为方程组的解。由以上两组计算可表明，该程序能满足追赶法解三对角线性方程《数值分析》数值实验XXX7组，其中要求系数矩阵满足要求，即非奇异，三对角且前（n-1）个顺序主子式都不为零。在又换了一组系数矩阵后，系数矩阵如下：A=[11000;22100;02310;00241;00025]即11221=23124125A，11=111d运行程序后出错，查找问题发现A进行Crout分解时，l1=1m2=2u1=1l2=a2-m2*u1=2-2×1=0再计算u2=c2/l2时无法计算，分析发现A的2阶顺序主子式为零，无法进行Crout分解，所以该程序还无法判断输入的系数矩阵的（n-1）阶顺序主子式是否为零，这是不足的地方，也无法判断系数矩阵是否是三对角矩阵，还需要在程序上进行修改完善。六实验总结1112212(1)1u(n1)()()1lumluALUlnmnln……………《数值分析》数值实验XXX8通过本次的课程实验设计，对追赶法求解三对角线性方程组有了深刻的认识和了解，用matlab编程实现了该解法；同时在实验过程中明白了系数矩阵需要满足一定的要求才能使用该种追赶法方法。