直接法、迭代法解线性方程组

1题目:用直接法、迭代法解线性方程组学生姓名：崔敬轩学号：2015210458专业班级：电气工程7班指导教师：王承竞2015年12月12日2一、直接法解线性方程组（选列主元高斯消元法）1高斯消元法具体地，解线性方程组AXb，先对系数矩阵A进行LU分解，即ALU，则AXbLUXb。令UXZ，则AXbLZb。利用1L解1ZLb，再通过1U解1XUZ。该方法讲A、b同时一起作初等行变换，优点：快捷、稳定；缺点：由硬件决定，运算规模受到限制。2选主元对线性方程组12011112xx它的解应为1211xx对于高斯消元法，0由特别小的正数替代，即1211112xx解得121112211(1)?xxx可以看到112(1)xx大小不能确定。经初等行变换使11a变成一个不接近0的较大的数，具体地，对于线性方程组12011112xx交换两方程的位置，即31211211xx解得2112121121xxxx3基于Matlab的数值试验3.1计算平台：①CPU（2.50GHz）、②RAM（4.00GB）、③OS（Windows7旗舰版）3.2参数设置：因为所用方法为直接法，所以不考虑算法误差、计算精度、迭代步数问题。3.3计算结果（程序代码详见附录Ⅰ）线性方程组1231231231132323110221xxxxxxxxx12123123226244385xxxxxxxx12342342342346224124221012821231426xxxxxxxxxxxxx方程组的解1231.00002.00001.0000xxx1230.62000.76000.0300xxx12341.00003.00002.00001.0000xxxx计算时间（s）0.0000000.0000000.000000⑴4⑵⑶53.4分析改进①对所求解的三个线性方程组，可以看到运算的时间都接近于0，说明用直接法解线性方程组的时间，尤其是低阶的线性方程组，所需时间很短。对于高阶、超高阶的线性方程组的运算能力有待进一步验证；②经Matlab计算的三个线性方程组的解，都为线性方程组的真实解。后面的尾数0为程序的默认保留小数位数；③改进：在选主元时，是以对应列中系数ija的绝对值最大的，作为主元的；如果先计算线性方程组每行的尺度（以每行中最大的||ija为尺度is），再通过ijias得比较来选择主元，效果会更好。6二、迭代法（Gauss-Seidel迭代法）解线性方程组1迭代法（分裂迭代）对线性方程组AXb，可找一个与A相近且易于求逆的矩阵Q，使原方程组等价于()QXQAXb11()XIQAXQb其中，I是单位阵。构造迭代式111()kkxIQAxQb若{kx}收敛，则*1*1()xIQAxQb分析{kx}的收敛性，1*1*1*||||||()()||||||||||kkkxxIQAxxIQAxx当且仅当1||||1IQA时，{kx}全局收敛。2Gauss-Seidel迭代当[]Q时，分裂迭代变为Gauss-Seidel迭代1111()kkkxDLxUxDb11111()()()[()]()kkkxDLUxDLbDLDLAxDLb3基于Matlab的数值试验3.1计算平台：①CPU（2.50GHz）、②RAM（4.00GB）、③OS（Windows7旗舰版）3.2参数设置：编写的m文件函数值误差限3110、最大迭代步数100M。3.3计算结果（程序代码详见附录Ⅱ）7线性方程组1231231231132323110221xxxxxxxxx12123123226244385xxxxxxxx12342342342346224124221012821231426xxxxxxxxxxxxx方程组的解1230.99831.99670.9976xxx1230.62000.76050.0298xxx12340.99542.99021.98591.0044xxxx迭代次数27776计算时间0.0468000.0312000.046800⑴⑵8⑶93.4分析改进①对于所求解的三个低阶方程组，求解时间相对直接法要长一些，但即便长一些，对于人的放映时间仍然是可以忽略的，随着运算阶数的增加，因为直接法硬件的限制及迭代法对内存要求不多的优点，猜测直接法的优势逐渐消失，而迭代法的应用价值会愈加明显；②线性方程组的解都接近于真实解，而不是真实解，这是所选误差限0.0001的缘故，如果进一步缩小误差限，会得到更加准确的解，并可能达到真实解，但现实生活中的运算大都出于工程实际需要，满足一定的精度要求即可；③改进：增大误差限，增大最大迭代步数M,求取所给线性方程组更准确的解。④无论是选列主元Gauss消元法，还是Gauss-Seidel迭代法，在举例计算中都只考虑了系数矩阵A的秩与增广矩阵（A,b）的秩相等且等于矩阵的长度n（即()(,)RARAbn）的情况，虽然在选列主元Gauss消元法中由()RA、(,)RAb、n三者间的关系对线性方程组的情况进行了判断，但对于()(,)RARAbn的情况，没有对线性方程组解的通式进行求解。10参考文献[1]同济大学数学系，线性代数（第五版）[M]，高等教育出版社，2007，P71-P78[2]DavidKincaid,WardCheney,数值分析（第三版）[M]，机械工业出版社，2005，P132-P136，P170-P172[3]徐跃良，数值分析[M],西南交通大学出版社，2005，P99-P103，P143-P14711附录Ⅰ选（列）主元Gauss消元法clear;clc;A=input('输入系数矩阵A:');b=input('输入常数项向量b(按行向量):');t0=cputime;%记录开始运算时间t0B=[Ab'];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;ifzhica0,disp('请注意：因为RA~=RB,所以此方程组无解！\n')returnendifRA==RBifRA==nfprintf('请注意：因为RA=RB=%d，所以此方程有唯一解！\n',n)X=zeros(n,1);forp=1:n-1t=find(abs(B(p:end,p))==max(abs(B(p:end,p))))+p-1;ifabs(B(t,p))~=abs(B(p,p))l=B(t,:);B(t,:)=B(p,:);B(p,:)=l;end%列主元判断fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endend%把系数矩阵A化为同解的上三角矩阵b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);end%从xn至x1逐个求解上三角方程组elsedisp('请注意：因为RA=RBn，所以此方程组有无穷多解！')returnendendt=cputime-t0;%计算程序运算时间tdisp('方程组的解为：');Xfprintf('\n运算时间t=%fs\n',t)12附录ⅡGauss-Seidel迭代法clear;clc;A=input('输入系数矩阵A=');b=input('输入常数项向量b=');t0=cputime;%记录开始运算时间t0N=length(b);%解向量的维数fprintf('库函数计算结果：');x=inv(A)*b%库函数计算结果x=zeros(N,1);%迭代初始值%-----(A=D-E-F)------D=diag(diag(A));E=-tril(A,-1);%下三角F=-triu(A,1);%上三角B=inv(D-E)*F;g=inv(D-E)*b;eps=0.0001;%相邻解的距离小于该数时，结束迭代%--------开始迭代-------fork=1:100%最大迭代次数为100fprintf('第%d次迭代：',k);y=B*x+g;fprintf('\n与上次计算结果的距离(2范数):%f\n',norm(x-y)^2);ifnorm(x-y)epsbreak;endx=yendt=cputime-t0;%计算程序运算时间txfprintf('\n运算时间t=%fs\n',t)