end
end
解方程组
(1)
在MATLAB工作窗口输入
>>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b)
运行后输出结果为
请注意:
因为RA=RB=n,所以此方程组有唯一解.
D=-0.1225
RA=3RB=3n=3
X=397.8654
-157.6242
-123.1120
解方程组
(2)
在MATLAB工作窗口输入
>>A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b)
运行后输出结果为
请注意:
因为RA=RB=n,所以此方程组有唯一解.
D=-762.0000
RA=4RB=4n=4
X=0.0000
-1.0000
1.0000
1.0000
LU分解法及MATLAB主程序为
functionhl=zhjLU(A)
[nn]=size(A);RA=rank(A);
D=det(A)
ifRA~=n
disp('请注意:
因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下:
'),RA,hl=det(A);
return
end
ifRA==n
forp=1:
n
h(p)=det(A(1:
p,1:
p));
end
hl=h(1:
n);
fori=1:
n
ifh(1,i)==0
disp('请注意:
因为A的r阶主子式等于零,所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:
'),hl;RA
return
end
end
ifh(1,i)~=0
disp('请注意:
因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:
')
forj=1:
n
U(1,j)=A(1,j);
end
fork=2:
n
fori=2:
n
forj=2:
n
L(1,1)=1;L(i,i)=1;
ifi>j
L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);
L(i,k)=(A(i,k)-L(i,1:
k-1)*U(1:
k-1,k))/U(k,k);
else
U(k,j)=A(k,j)-L(k,1:
k-1)*U(1:
k-1,j);
end
end
end
end
hl;RA,U,L
end
end
解方程组
(1)
在MATLAB工作窗口输入
>>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];h1=zhjLU(A)
运行输出结果为
请注意:
因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:
D=9.8547
RA=3
U=3.01006.03001.9990
04.1600-2.0734
005.3016
L=1.000000
0.42191.00000
0.3279-1.63161.0000
h1=3.01004.8635-0.1225
解方程组
(2)
在MATLAB工作窗口输入
>>A=[10-701;-32.09999962;5-15-1;2102];h1=zhjLU(A)
运行后输出结果为
请注意:
因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:
D=-762.0000
RA=4
U=10.0000-7.000001.0000
02.10006.00002.3000
00-2.1429-4.2381
0-0.0000012.7333
L=1.0000000
-0.30001.000000
0.50001.19051.0000-0.0000
0.20001.14293.20001.0000
h1=10.0000-0.0000-150.0001-762.0001
(2)在MATLAB工作窗口输入
>>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];A(1,1)=3;A(1,3)=0.990;[RA,RB,n,X]=liezhu(A,b)
请注意:
因为RA=RB=n,所以此方程组有唯一解.
RA=3RB=3n=3
X=-4.0264
1.9193
1.5210
hi=3.00004.82199.8547
在MATLAB工作窗口输入
x=[397.8654;-157.6242;-123.1120]';x1=[-4.0264;1.9193;1.5210]';wucha=x1-x
运行后输出结果为
wucha=-401.8918159.5435124.6330
(3)在MATLAB工作窗口输入
>>A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;b(2,1)=5.9;b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b)
运行后输出结果为
请注意:
因为RA=RB=n,所以此方程组有唯一解.
RA=4RB=4n=4
X=0.0000
-1.0000
1.0000
1.0000
h1=10.0000-0.0000-150.0000-762.0000
在MATLAB工作窗口输入
>>x=[0;-1;1;1]';x1=[0;-1;1;1]';wucha=x1-x
运行后输出结果为
wucha=0000
(4)解方程组
(1)
在MATLAB工作窗口输入
>>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];B=inv(A)
运行后结果为
B=-268.9293538.3418128.4529
106.7599-213.4281-50.9561
83.3992-166.8022-39.7090
在MATLAB工作窗口输入
>>b=[1;1;1];x=inv(A)*b
运行后结果为
x=397.8654
-157.6242
-123.1120
在MATLAB工作窗口输入
>>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];A(1,1)=3;A(1,3)=0.990;B=inv(A)
运行输出结果为
B=3.3424-6.1983-1.1705
-1.32692.74420.5020
-1.03652.06820.4893
在MATLAB工作窗口输入
>>b=[1;1;1];x=inv(A)*b
运行后输出结果为
x=-4.0264
1.9193
1.5210
解方程组
(2)
在MATLAB工作窗口输入
>>A=[10-701;-32.09999962;5-15-1;2102];B=inv(A)
运行后结果为
B=-0.0223-0.09840.11810.1686
-0.1601-0.11810.14170.2690
0.01080.10630.0724-0.0755
0.10240.1575-0.18900.1969
在MATLAB工作窗口输入
>>b=[8;5.900001;5;1];x=inv(A)*b
运行后输出结果为
x=0
-1.0000
1.0000
1.0000
在MATLAB工作窗口输入
>>A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;B=inv(A)
运行后输出结果为
B=-0.0223-0.09840.11810.1686
-0.1601-0.11810.14170.2690
0.01080.10630.0724-0.0755
0.10240.1575-0.18900.1969
在MATLAB工作窗口输入
>>b=[8;5.900001;5;1];b(2,1)=5.9;x=inv(A)*b
运行后输出结果为
x=-0.0000
-1.0000
1.0000
1.0000
五、实验结果分析:
实验的数学原理很容易理解,也容易上手。
把运算的结果带入原方程组,可以发现符合的还是比较好。
这说明列主元消去法计算这类方程的有效性。
当A可逆时,能够将计算进行到底,列主元法就能确保算法的稳定,而且计算量不大。
直接三角消去过程,实质上是将A分解为两个三角矩阵的乘积A=LU,并求解Ly=b的过程。
回带过程就是求解上三角方程组Ux=y。
所以在实际的运算中,矩阵L和U可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法。
通过以上的计算比较,方程组
(1)具有严重的病态性。
当系数矩阵有微小的变化时,wucha=-401.8918159.5435124.6330,所得的解与原方程组的解有很大的相对误差。
方程组
(2)中当系数矩阵A和b有微小变化时,wucha=0000,所得的解与方程组的解没有相对误差。
所以方程组
(2)是良性的。
用MATLAB内部函数inv通过求逆矩阵,然后通过x=inv(A)*b也可以求出方程组的解,但是没有列主元高斯消去法具有良好的稳定性。
det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同,具有方便快捷的优点。