数值计算实验|解线性方程组的消元法

一、实验目的

本实验主要涉及解线性方程组的直接法中的消元法和列主元消去法。通过上机计算使学生学会程序的录入和Matlab的使用与操作；通过计算结果的比较，使学生了解求线性方程组的精度不但与方法有关，而且与问题的性态有关。

二、实验内容

（1）算法 （2）程序 ●程序代码

function x=shiyan1(A,b) [n,n]=size(A); for k=1:n p(k)=k; end for k=1:n-1 %Ñ¡ÁÐÖ÷Ôª %选列主元 [temp,maxk]=max(abs(A(k:n,k))); maxk=maxk+k-1; if(maxk~=k) temp1=A(k,1:n); A(k,1:n)=A(maxk,1:n); A(maxk,1:n)=temp1; temp2=p(k); p(k)=p(maxk); p(maxk)=temp2; end %¼ÆËãPA=LU %计算PA=LU if(A(k,k)~=0) A(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); end end clear temp temp1 temp2?%ÐγÉÓÒ¶ËÏî %形成右端项 b=b(p); %Çó½âLy=b %求解Ly=b for k=2:n b(k)=b(k)-A(k,1:k-1)*b(1:k-1); end %Çó½âUx=y %求解Ux=y b(n)=b(n)/A(n,n); for k=n-1:-1:1 b(k)=(b(k)-A(k,k+1:n)*b(k+1:n))/A(k,k); end %加入赋值x，保证x被赋值 for k=1:n x(k)=b(k); end

●程序截图

三、实验任务

（1）将上述程序录入计算机，并进行调试。 调试成功。 （2）用调试好的程序解决如下问题： 用Gauss列主元素消去法求Ax=b的解，其中 ●（a）题命令行窗口程序输入与输出结果如下：

>> A=[1 -1 2 -1;2 -2 3 -3;1 1 1 0;1 -1 4 3] A = 1 -1 2 -1 2 -2 3 -3 1 1 1 0 1 -1 4 3 >> b=[-8 -20 -2 4]' b = -8 -20 -2 4 >> format long; >> x=shiyan1(A,b) x = -6.999999999999995 2.999999999999998 1.999999999999997 2.000000000000002

●（a）题窗口截图如下： ●（b）题命令行窗口程序输入与输出结果如下：

>> A=[2.0000 -2.0000 3.0000 -3.0000;0.5000 2.0000 -0.5000 1.5000;0.5000 0 2.5000 4.5000;0.5000 0 0.2000 -0.4000] A = 2.000000000000000 -2.000000000000000 3.000000000000000 -3.000000000000000 0.500000000000000 2.000000000000000 -0.500000000000000 1.500000000000000 0.500000000000000 0 2.500000000000000 4.500000000000000 0.500000000000000 0 0.200000000000000 -0.400000000000000 >> b=[-14 6 4 4]' b = -14 6 4 4 >> format long; >> x=shiyan1(A,b) x = 25.161290322580641 -16.612903225806445 -22.129032258064509 10.387096774193546

●（b）题窗口截图如下： （3）调用Matlab中的“\”解上述算例。 ●（a）题命令行窗口程序输入与输出结果如下：

>> A=[1 -1 2 -1;2 -2 3 -3;1 1 1 0;1 -1 4 3] A = 1 -1 2 -1 2 -2 3 -3 1 1 1 0 1 -1 4 3 >> b=[-8 -20 -2 4]' b = -8 -20 -2 4 >> format long; >> x=A\b x = -6.999999999999998 3.000000000000000 1.999999999999999 2.000000000000000 >> format long; >> x=shiyan1(A,b) x = -6.999999999999995 2.999999999999998 1.999999999999997

●（a）题窗口截图如下：

●（b）题命令行窗口程序输入与输出结果如下：

>> A=[2.0000 -2.0000 3.0000 -3.0000;0.5000 2.0000 -0.5000 1.5000;0.5000 0 2.5000 4.5000;0.5000 0 0.2000 -0.4000] A = 2.000000000000000 -2.000000000000000 3.000000000000000 -3.000000000000000 0.500000000000000 2.000000000000000 -0.500000000000000 1.500000000000000 0.500000000000000 0 2.500000000000000 4.500000000000000 0.500000000000000 0 0.200000000000000 -0.400000000000000 >> b=[-14 6 4 4]' b = -14 6 4 4 >> format long; >> x=A\b x = 25.161290322580644 -16.612903225806452 -22.129032258064516 10.387096774193548 >> format long; >> x=shiyan1(A,b) x = 25.161290322580641 -16.612903225806445 -22.129032258064509 10.387096774193546

●（b）题窗口截图如下：

（4）分别用上述程序和Matlab中的“\”求线性方程组 ①当n=10时： ●用程序解：

>> A=hilb(10); b=(1:10)'; format long; shiyan1(A,b)' ans = 1.0e+08 * -0.000009998189333 0.000979943729503 -0.023281474448041 0.233002589296527 -1.210665589642735 3.594195329722782 -6.322467936951993 6.510563619029391 -3.622833473596718 0.840567489799199

●窗口截图如下：

●用“\”解：

>> A=hilb(10); b=(1:10)'; format long; A\b ans = 1.0e+08 * -0.000009998316280 0.000979953705980 -0.023281671354817 0.233004270001089 -1.210673190038930 3.594215286826776 -6.322499397976280 6.510592970934195 -3.622848408406987 0.840570683406235

●窗口截图如下：

②当n=20时： ●用程序解：

>> A=hilb(20); b=(1:20)'; format long; shiyan1(A,b)' ans = 1.0e+11 * -0.000000019403479 0.000002224327300 -0.000052948754714 0.000292943938818 0.003055979438376 -0.048410770059000 0.272233814835646 -0.801741935096167 1.279514831483144 -1.002230606371555 0.512286794841114 -1.217355338178257 1.628287819432132 1.072830291206492 -4.445149807227658 4.240040585069163 -2.062975568195349 1.022957695111095 -0.640594667648723 0.187009007403822

●窗口截图如下：

●用“\”解：

>> A=hilb(20); b=(1:20)'; format long; A\b 警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND = 5.231543e-20。 ans = 1.0e+12 * 0.000000026996781 -0.000004291664839 0.000167935154054 -0.002812084432282 0.024839697941581 -0.127222322914451 0.389583742850322 -0.691767340721250 0.624587527253181 -0.263848481338642 0.664265987209912 -1.919935436473158 2.166604558474094 -1.194549954842270 1.119552208472303 -1.423463283478544 0.113473492323463 1.343704151687982 -1.094214762919976 0.271038680397702

●窗口截图如下： ③当n=30时： ●用程序解：

>> A=hilb(30); b=(1:30)'; format long; shiyan1(A,b)' ans = 1.0e+11 * -0.000000076924758 0.000010845649352 -0.000357391268411 0.004551869910115 -0.023557638195059 0.003581210206813 0.519662034582864 -2.541090293636888 5.709794455735597 -6.689213142862793 4.969593290464653 -6.166956872577877 7.479141272547461 -0.018684245719832 -4.649734224423247 -2.515216624415138 4.017514943934164 -1.998991797957529 7.340894914261927 -4.295019873949639 -0.761317461336779 -4.074440186140318 2.577440068806253 -1.009537797519588 8.047431197540366 -8.960421097157296 5.484779935732302 -4.493095803645594 2.479166628311472 -0.435927375365802

●窗口截图如下：

●用“\”解：

>> A=hilb(30); >> b=(1:30)'; format long; A\b 警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND = 1.647474e-19。 ans = 1.0e+11 * -0.000000046948880 0.000007445182443 -0.000286845711456 0.004647140967376 -0.038725217332101 0.179537454343672 -0.457395748443258 0.516862593815496 0.179335577967648 -0.936842465077658 0.123278466926140 1.519340914315617 -1.846475806058420 1.932670121041188 -3.158502173614151 3.357663662803517 -1.066177187175942 -0.774506960047759 -1.006843372806605 1.789951172911067 1.220893536364581 -0.499394753292967 -1.525998990848227 -0.346857476965359 -1.376817806491880 3.020536079430076 0.914669212972973 -1.628824403751380 -0.771123067991843 0.675379681756272

●窗口截图如下：

（5）对上述算例的计算结果进行分析，你有哪些启示？ 对于第（4）小问中的，当n=20和n=30时：两种方法求得的值不相同，用“\”求解出来的结果严重失真，甚至在MATLAB程序里出现了警告，而用Gauss求得准确解。因为Hilbert矩阵是“病态矩阵”，所以会导致所得的解不准确。 所以我们在使用程序时，一定要注意警告里的提示，虽然它并不能阻止程序的运行，但是它可能反映了问题与结果的准确性，虽然程序可以很快捷地求出解，但是如果我们没有掌握好理论知识而疏忽程序的编写的话，可能反而会导致事倍功半。 此外，程序的编写也不是只有唯一解答的，要不断掌握新的知识，这样才能尽可能用简单的方法求出复杂的问题。