MATLAB实验一解线性方程组的直接法.docx
《MATLAB实验一解线性方程组的直接法.docx》由会员分享，可在线阅读，更多相关《MATLAB实验一解线性方程组的直接法.docx（20页珍藏版）》请在冰豆网上搜索。
MATLAB实验一解线性方程组的直接法
实验报告
课程名称数值分析
实验项目解线性方程组的直接法
专业班级姓名学号
指导教师成绩日期月曰
1.实验目的
1掌握程序的录入和matlab的使用和操作；
2、了解影响线性方程组解的精度的因素——方法与问题的性态。
3、学会Matlab提供的“”的求解线性方程组。
2.实验要求
1按照题目要求完成实验内容；
2、写出相应的Matlab程序；
3、给出实验结果（可以用表格展示实验结果）；
4、分析和讨论实验结果并提出可能的优化实验。
5、写出实验报告。
3.实验步骤
1用LU分解及列主元高斯消去法解线性方程组
z10
-7
0
1、
1
f、
X1
，z8'
-3
2.099999
6
21
X2
5.900001
a）
—
5
-1
5
一1
X3
5
<2
1
0
2丿k丿
<1丿
输出Ax二b中系数A二LU分解的矩阵L和U，解向量x和det（A）;用列主元法的行交换次序解向量x和求det（A）;比较两种方法所得结果。
2、用列主高斯消元法解线性方程组Ax二b
p.01
6.03
1.99、
p1、
I
（1）、
1.27
4.16
-1.23
x
=
1
0987
-4.81
9.34‘
Z丿
l1丿
广3.00
6.03
1.99、
■行
（2）、
1.27
4.16
-1.23
X2
=
1
0990
-4.81
9.34”
W丿
J丿
分别输出A,b,det（A）,
解向量x,
（1）
中A的条件数。
分析比较
（1）、
（2）的
计算结果
cond（A）2.若令
求解（A、：
A）（x，x）二b，输出向量x和、胡2，从理论结果和实际计算两方面分析线
分量的有效位数如何随n变化，它与条件数有何关系？
当n多大时x连一位有效数字也
3、线性方程组Ax=b的A和b分别为
n=6时：
没有了？
GAUSS列主消去法求得的x
x的有效数字
将每种情形的两个结果进行表格对比，如:
四、实验结果
五、讨论分析
（对上述算例的计算结果进行比较分析，主要说清matlab的算符与消去法的适
用范围不同，自己补充）
六、改进实验建议
（自己补充）
1•列主元的高斯消去法
利用列主元的高斯消去法matlab程序源代码：
首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。
functionx=gaussMethod（A,b）
%高斯列主元消去法，要求系数矩阵非奇异的，%
n=size（A,1）;
ifabs（det（A））<=1e-8
error（'系数矩阵是奇异的'）;
return;
end
%
fork=1:
n
ak=max（abs（A（k:
n,k）））;
index=find（A（:
k）==ak）;
iflength（index）==0
index=find（A（:
k）==-ak）;
end
%交换列主元
temp=A（index,:
）;
A（index,:
）=A（k,:
）;
A（k,:
）=temp;
temp=b（index）;b（index）=b（k）;b（k）=temp;
%消元过程
fori=k+1:
nm=A（i,k）/A（k,k）;
%消除列元素
A（i,k+1:
n）=A（i,k+1:
n）-m*A（k,k+1:
n）;b（i）=b（i）-m*b（k）;
end
end
%回代过程
x（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1;
x（k）=（b（k）-A（k,k+1:
n）*x（k+1:
n）'）/A（k,k）;
end
x=x';
end
然后调用gaussMethod函数，来实现列主元的高斯消去法。
在命令框中输入下列命令：
CommandWindow
»clear
»A=[10,-7,031；-X2.099999,6,2;5,-1,5,-1;2,1,0,2];
b=[8.5,900001；5；1]；
k=gaussMethod（A3b）
detA=det（A）
输出结果如下：
0.0000-LOOOO1.0000LOOOO
detA=
^762*0001
利用LU分解法及matlab程序源代码：
function[L,U]=myLU（A）%实现对矩阵A的LU分解，L为下三角矩阵A[n,n]=size（A）;
L=zeros（n,n）;
U=zeros（n,n）;
fori=1:
n
L（i,i）=1；
end
fork=1:
n
forj=k:
n
U（k,j）=A（k,j）-sum（L（k,1:
k-1）.*U（1:
k-1,j）'）;end
fori=k+1:
nL（i,k）=（A（i,k）-sum（L（i,1:
k-1）.*U（1:
k-1,k）'））/U（k,k）;end
end
在命令框中输入下列命令：
>>a=[10-701;-32.09999962;5-15-1;2102]a=
10.0000
-7.0000
0
1.0000
-3.0000
2.1000
6.0000
2.0000
5.0000
-1.0000
5.0000
-1.0000
2.0000
1.0000
0
2.0000
>>
l=
[l,u]=lu（a）
1
1.0000
0
0
0
-0.3000
-0.0000
1.0000
0
0.5000
1.0000
0
0
0.2000
0.9600
-0.8000
1.0000
u=
10.0000
-7.0000
0
1.0000
0
2.5000
5.0000
-1.5000
0
0
6.0000
2.3000
0
0
0
5.0800
>>
b=[85.90000151]'
b=
8.0000
5.9000
5.0000
1.0000>>y=l\by=
8.0000
1.0000
8.3000
5.0800
>>x1=U\xx1=
0.0000
-1.0000
1.0000
1.0000
>>det仁det（a）
detl=-762.0001
2、
（1）在MATLAB窗口：
>>A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34]A=
3.0100
6.0300
1.9900
1.2700
4.1600
-1.2300
0.9870
-4.8100
9.3400
>>b=[111]'
b=
1
1
1
>>[x1,det1,index]=Gauss（A,b）x1=
1.0e+03*
1592.599624841381
-631.9113762025488
-493.6177247593899
det1=
-0.0305
index=
1
⑵在MATLAB窗口：
>>A=[3.006.031.99;1.274.16-1.23;0.990-4.819.34]A=
3.0000
6.0300
1.9900
1.2700
4.1600
-1.2300
0.9900
-4.8100
9.3400
>>b=[111]'
b=
1
1
1
>>[x2,det2,index]=Gauss5555（A,b）
x2=
119.5273
-47.1426
-36.8403
det2=
-0.4070index=
1
3、在MATLAB窗口：
A=[10787;7565;86109;75910];
b=[32233331]';
x=A\b
b仁[32.122.933.130.9]';
x1=A\b1
A1=[1078.17.2;7.085.0465;85.989.899;6.99599.98];x2=A1\b
delta_b=norm（b-b1）/norm（b）
delta_A=norm（A-A1）/norm（A）
delta_x1=norm（x-x1）/norm（x）
deltax2=norm（x-x2）/norm（x）
cond_A=cond（A）
x=
1.0000
1.0000
1.0000
1.0000x1=
9.2000-12.6000
4.5000
-1.1000x2=
-9.5863
18.3741-3.2258
3.5240delta_b=
0.0033delta_A=
0.0076delta_x1=
8.1985delta_x2=10.4661cond_A=
2.9841e+03
3、在MATLAB窗口：
A=[10787;7565;86109;75910];
b=[32233331]';
x=A\b
b1=[32.122.933.130.9]';
x1=A\b1
A1=[1078.17.2;7.085.0465;85.989.899;6.99599.98];x2=A1\b
deltab=norm（b-b1）/norm（b）
delta_A=norm（A-A1）/norm（A）delta_x1=norm（x-x1）/norm（x）delta_x2=norm（x-x2）/norm（x）cond_A=cond（A）
x=
1.0000
1.0000
1.0000
1.0000
x1=
9.2000
-12.6000
4.5000
-1.1000
x2=
-9.5863
18.3741-3.2258
3.5240
delta_b=
0.0033
delta_A=
0.0076
delta_x1=
8.1985
delta_x2=
10.4661cond_A=
2.9841e+03
4、
k=13
forn=2:
6
a=hilb（n）;
co（n）=cond（a,inf）;endx=1:
6;
plot（x,co）;
b=zeros（k）;
x1仁b;
x0=b;
r=b;
fori=2:
k
x=（linspace（1,1,i））';
x0（1:
i,（i-1））=x;
H=hilb（i）;
b0=H*x;
b（1:
i,（i-1））=b0;
x1=gauss2（H,bO）;
r（1:
i,（i-1））=b0-H*x1;
x11（1:
i,（i-1））=x1;
end
dx=x11-x0;
结果如下：
co=[0,27,748,28375,943656,29070279]可见，条件数随着n的增大，急剧增加，如图1所
示。
将
（2）求得的结果（dx,x11）整理得
n=2时：
x
△x
rn
有效数字
1
4.44E-16
0
16
1
-6.66E-16
0
15
n=3时:
x
△x
rn
有效数字
1
-1.33E-15
-2.22E-16
15
1
9.55E-15
0
14
1
-9.99E-15
0
14
n=4时:
x
△x
rn
有效数字
1
-2.35E-14
0
14
1
2.56E-13
0
13
1
-6.11E-13
1.11E-16
12
1
3.96E-13
0
13
n=5时:
x
△x
rn
有效数字
1
-1.21E-14
4.44E-16
14
1
6.97E-14
0
13
1
1.18E-13
2.22E-16
13
1
-6.17E-13
0
12
1
4.59E-13
-1.11E-16
13
n=6时:
x
△x
rn
有效数字
1
-9.28E-13
0
12
1
2.67E-11
0
11
1
-1.82E-10
0
10
1
4.75E-10
0
10
1
-5.26E-10
0
9
1
2.08E-10
0
10
n=7时:
x
△x
rn
有效数字
1
-9.26E-12
0
11
1
3.71E-10
0
10
1
-3.59E-09
2.22E-16
9
1.00000001
1.40E-08
0
8
0.99999997
-2.57E-08
0
8
1.00000002
2.22E-08
1.11E-16
8
0.99999999
-7.30E-09
0
8
n=8时:
x
△x
rn
有效数字
1
-2.82E-11
4.44E-16
11
1
1.53E-09
0
9
0.99999998
-2.01E-08
2.22E-16
8
1.00000011
1.09E-07
0
7
0.9999997
-2.95E-07
-2.22E-16
7
1.00000042
4.19E-07
-1.11E-16
7
0.9999997
-2.99E-07
1.11E-16
7
1.00000008
8.44E-08
0
7
n=9时:
x
△x
rn
有效数字
1
-2.75E-10
4.44E-16
10
1.00000002
1.90E-08
0
8
0.99999968
-3.20E-07
0
7
1.00000228
2.28E-06
4.44E-16
6
0.99999164
-8.36E-06
0
5
1.00001706
1.71E-05
0
5
0.99998042
-1.96E-05
-1.11E-16
5
1.00001182
1.18E-05
0
5
0.99999708
-2.92E-06
0
6
n=10时:
x
△x
rn
有效数字
1
-9.05E-10
4.44E-16
9
1.00000008
7.76E-08
0
7
0.99999835
-1.65E-06
0
6
1.00001491
1.49E-05
4.44E-16
5
0.99992911
-7.09E-05
0
4
1.00019443
0.000194426
0
4
0.99968164
-0.000318365
-1.11E-16
4
1.00030715
0.000307149
0
4
0.99983898
-0.000161019
0
4
1.00003537
3.54E-05
0
5
n=11时:
x
△x
rn
有效数字
0.999999995
-5.16E-09
0
8
1.000000539
5.39E-07
-4.44E-16
6
0.999986041
-1.40E-05
0
5
1.000155875
0.00015588
0
4
0.999072325
-0.0009277
0
3
1.003258718
0.00325872
0
3
0.992910161
-0.0070898
0
2
1.009658807
0.00965881
-2.22E-16
2
0.991981908
-0.0080181
0
3
1.003707654
0.00370765
0
3
0.999267974
-0.000732
-2.22E-16
3
n=12时:
x
△x
rn
有效数字
0.999999965
-3.49E-08
8.88E-16
8
1.000004409
4.41E-06
4.44E-16
6
0.999861744
-0.0001383
0
4
1.001879711
0.00187971
0
3
0.986241983
-0.013758
0
2
1.060375974
0.06037597
2.22E-16
2
0.831939173
-0.1680608
2.22E-16
1
1.303973363
0.30397336
0
1
0.643862006
-0.356138
1.11E-16
1
1.260684018
0.26068402
2.22E-16
1
0.891666924
-0.1083331
1.11E-16
1
1.019510753
0.01951075
0
2
n=13时:
x
△x
rn
有效数字
1.000000069
6.89E-08
8.88E-16
7
0.999989224
-1.08E-05
-4.44E-16
5
1.000413538
0.00041354
-2.22E-16
4
0.993147495
-0.0068525
2.22E-16
3
1.061246208
0.06124621
-2.22E-16
2
0.669261244
-0.3307388
2.22E-16
1
2.149074131
1.14907413
0
0
-1.65417119
-2.6541712
0
0
5.11854808
4.11854808
-1.11E-16
0
-3.243049939
-4.2430499
1.11E-16
0
3.783004896
2.7830049
0
0
-0.051792817
-1.0517928
1.11E-16
0
1.174329117
0.17432912
1.11E-16
1
五、分析讨论:
实验的数学原理很容易理解，也容易上手。
把运算的结果带入原方程组，可以发现符合的还是比较好。
这说明列主元消去法计算这类方程的有效性。
当A可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。
直接三角消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU并求解Ly=b的过程。
回带过程就是求解上三角方程组Ux=yo所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。
通过以上的计算比较，2.题方程组具有严重的病态性。
当系数矩阵有微小的变化时，wucha=-401.8918159.5435124.6330，所得的解与原方程组的解有很大的相对误差。
1题方程组中当系数矩阵A和b有微小变化时，wucha=0000，所得的解与方程组的解没有相对误差。
所以1题方程组是良性的。
用MATLA内部函数inv通过求逆矩阵，然后通过x=inv（A）*b也可以求出方程组的解，但是没有列主元高斯消去法具有良好的稳定性。
det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有方便快捷的优点。
题四可以看出，条件数越大，有效位数越少，当n=13时，出现有效位数为0的情况。
附：
高斯列主消去法源程序代码
function[x,det,index]=Gauss（A,b）
%?
o?
?
D?
•?
3iXe^?
aD?
-?
aGauss?
?
OY•…£
?
?
?
?
D
%A---
方程组矩阵
%b---
方程组右端
%x---
方程组的解
%det---
方程组行列式
%index---
index=0表示求解失败，index=1
表示求解成功
[n,m]=size（A）;nb=length（b）;ifn~=m
error（'TherowsandcolumnsofmatrixAmustbeequal!
'）;return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthedimensionofb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
%选主元
a_max=0;
fori=k:
n
ifabs（A（i,k））>a_max
a_max=abs（A（i,k））;r=i;
end
end
ifa_max<1e-20
index=0;
return;
end
%交换两行
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
%消元过程
fori=k+1:
nm=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
endb（i）=b（i）-m*b（k）;enddet=det*A（k,k）;
end
det=det*A（n,n）;
%回代过程
ifabs（A（n,n））<1e-20
index=0;
return;
end
fork=n:
-1:
1
forj=k+1:
nb（k）=b（k）-A（k,j）*x（j）;endx（k）=b（k）/A（k,k）;
end
升级会员 