一、实验目的
本实验主要涉及解线性方程组的直接法中的消元法和列主元消去法。通过上机计算使学生学会程序的录入和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矩阵是“病态矩阵”,所以会导致所得的解不准确。
所以我们在使用程序时,一定要注意警告里的提示,虽然它并不能阻止程序的运行,但是它可能反映了问题与结果的准确性,虽然程序可以很快捷地求出解,但是如果我们没有掌握好理论知识而疏忽程序的编写的话,可能反而会导致事倍功半。
此外,程序的编写也不是只有唯一解答的,要不断掌握新的知识,这样才能尽可能用简单的方法求出复杂的问题。