数值分析线性方程组的实验报告包含代码解析Word文件下载.docx
《数值分析线性方程组的实验报告包含代码解析Word文件下载.docx》由会员分享,可在线阅读,更多相关《数值分析线性方程组的实验报告包含代码解析Word文件下载.docx(18页珍藏版)》请在冰豆网上搜索。
(i=1,2,3……,n).将(2-1)中第i个方程减去第一个方程乘以
(i=1,2,3……,n),完成第一次消元,
(2-2)
其中:
;
:
简记为:
按上述方法完成n-1次消元后得到。
同解的三角方程组:
2回代过程:
按逆序逐步回代得到方程的解。
(3)算法:
(5)Matlab程序代码:
function[RA,RB,n,X]=gaus(A,b)
B=[A,b];
n=length(b);
RA=rank(A);
RB=rank(B);
zhicha=RB-RA;
if(zhicha~=0)
disp('
因为RA≠RB,所以此方程组无解'
);
return;
end
ifRA==RB
ifRA==n
因为RA=RB=n,所以此方程有唯一解'
X=zeros(n,1);
forp=1:
n-1
fork=p+1:
n
m=B(k,p)/B(p,p);
B(k,p:
n+1)=B(k,p:
n+1)-m*B(p,p:
n+1);
end
b=B(1:
n,n+1);
A=B(1:
n,1:
n);
X(n)=b(n)/A(n,n);
forq=n-1:
-1:
1
X(q)=(b(q)-sum(A(q,q+1:
n)*X(q+1:
n)))/A(q,q);
else
因为RA=RB<
n,所以此方程有无穷解'
运行检测:
>
A=[111;
12-33;
-183-1];
b=[6;
15;
-15];
[RA,RB,n,X]=gaus(A,b)
因为RA=RB=n,所以此方程有唯一解
RA=
3
RB=
n=
X=
1.0000
2.0000
3.0000
二、主元素法
1.列主元素法消元
(1)基本思想:
在每次消元前,在要消去未知数的系数中找到绝对值最大的系数作为主元,通过对换行将其换到对角线上,然后进行消元.
(2)消元过程:
与Gauss很类似,每次对对角线换成最大的值,后面过程与Gauss基本相同。
如此经过n-1步,(2-1)的增广矩阵[A,b]被化为上三角矩阵;
回代过程:
同Gauss算法一样回代求解。
det<
-1
对于k=1,2,3,…n-1;
|aik,k|=max|aik|(k≤i≤n)
(i)如果aik=0,则停止计算(det(A)=0)
(ii)如果aik=k,则zhuan(i)
换行akj<
-->
aik,j(j=k,k+1,…n)
bk<
bik
后边就是Gauss消元
(4).Matlab程序
functionX=liezhu(A,b)
n=length(b);
RA=rank(A);
RB=rank(B);
ifRA==n
因为系数矩阵与增广矩阵的秩相同且为n,此方程有唯一解'
)
forp=1:
n-1
[y,j]=max(abs(B(p:
n,p)));
C=B(p,:
B(p,:
)=B(j+p-1,:
B(j+p-1,:
)=C;
fork=p+1:
m=B(k,p)/B(p,p);
end
b=B(1:
A=B(1:
X(n)=b(n)/A(n,n);
forq=n-1:
else
因为系数矩阵与增广矩阵的秩相同且小于n,方程组有无穷解'
)
else
因为系数矩阵与增广矩阵的秩不相同,所以此方程组无解'
测试结果:
clearall
clc
b=[615-15]'
;
liezhu(A,b)
因为系数矩阵与增广矩阵的秩相同且为n,此方程有唯一解
ans=
3.0000
例2测试结果
clearall
A=[0.501.13.1;
2.04.50.36;
5.00.966.5];
b=[6.00.0200.96]'
-2.6000
2.0000
三、直接三角分解法
1.高斯消元的矩阵表示
第一次消元
经过K次消元后得到:
则
(1)LU分解:
,
因为A的各元素已知,可求出L和U的各元素。
(2)看A的第k行,有
,再看A的第k列,有
。
这样就可以计算计算出U的第k行和L的第k列的元素。
(3)从k=1到k=n交替计算U和L,就能逐步计算出U和L的全部元素。
(4)分两步求解LUx=b,先解方程组Ly=b,因为L是下三角矩阵,求解只要逐步向前回代。
第二步是解方程组Ux=y,因为U是上三角矩阵,用逐次向后回代的方法求出x。
3.Matlab程序
function[L,U,Y,X]=Doolittle(A,b)
L=zeros(n,n);
U=zeros(n,n);
%先将A进行LU分解
fori=1:
%把第一列第一行求出
forj=i:
ifi==1
A(i,j)=A(i,j);
ifj>
i
A(j,i)=A(j,i)/A(i,i);
%下面计算U的第i行值
ifi~=1
m=0;
fork=1:
i-1
m=m+A(i,k)*A(k,j);
A(i,j)=A(i,j)-m;
%下面计算L的第i列值
ifi~=1&
&
j>
l=0;
l=l+A(j,k)*A(k,i);
A(j,i)=(A(j,i)-l)/A(i,i);
%将LU的值分别加入
U(i,j)=A(i,j);
forj=1:
ifi==j
L(i,j)=1;
L(i,j)=A(i,j);
%LY=B,UX=Y,进行求解
Ab=[L,b]
Y=zeros(n,1);
Y
(1)=Ab(1,n+1)/Ab(1,1);
fori=2:
Y(i)=Y(i)+Ab(i,j)*Y(j);
Y(i)=(Ab(i,n+1)-Y(i))/Ab(i,i);
Ab=[U,Y];
X=zeros(n,1);
X(n)=Ab(n,n+1)/Ab(n,n);
fori=n-1:
forj=i+1:
X(i)=X(i)+Ab(i,j)*X(j);
X(i)=(Ab(i,n+1)-X(i))/Ab(i,i);
[L,U,Y,X]=Doolittle(A,b)
运行结果:
Ab=
1.0000006.0000
12.00001.0000015.0000
-18.0000-1.40001.0000-15.0000
L=
1.000000
12.00001.00000
-18.0000-1.40001.0000
U=
1.00001.00001.0000
0-15.0000-9.0000
004.4000
Y=
6.0000
-57.0000
13.2000
四、解三角方程组的追赶法
1.设系数矩阵为三对角矩阵
则方程组Ax=f称为三对角方程组。
2.设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记
先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。
求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法吏为紧凑。
其计算公式为:
3.算法:
1.输入a=(a2,…,an),b=(b1,…,bn),c=(c1,…,cn-1),d=(d1,…,dn),n
2.对i=2,3,…,n
ai/bi-1=>
ai(li)
bi-ci-1ai=>
bi(ui)
di-aidi-1=>
di(yi)
3.dn/bn=>
dn(xn)
4.对i=n-1,n-2…,1
(di-cidi+1)/bi=>
di(xi)
5.输出d=x,停机
4.Matlab程序:
clc;
clear;
%如果需要验证其他矩阵的解,只需更改A、B和n的值即可
A=[2,1,0,0;
1,3,1,0;
0,1,4,-1;
0,0,2,3]
B=[3;
5;
4;
5]
n=4;
%获得A中不为0的数存到相应数组中
b(i)=A(i,j);
ifi==j-1
c(j-1)=A(i,j);
ifj==i-1
a(i)=A(i,j);
%将变换后的值存入相应矩阵中
u
(1)=b
(1);
n;
l(i)=a(i)/u(i-1);
u(i)=b(i)-l(i)*c(i-1);
U(i,j)=u(i);
U(i,j)=c(i);
L(i,j)=l(i);
disp('
增广矩阵'
AB=[L,B]
B=zeros(n,1);
B
(1)=AB(1,n+1)/AB(1,1);
B(i)=B(i)+AB(i,j)*B(j);
B(i)=(AB(i,n+1)-B(i))/AB(i,i);
解方程LY=B得'
Y=B
AB=[U,B];
B(n)=AB(n,n+1)/AB(n,n);
方程的解为:
'
X=B
A=
2100
1310
014-1
0023
B=
5
4
增广矩阵
AB=
1.00000003.0000
0.50001.0000005.0000
00.40001.000004.0000
000.55561.00005.0000
解方程LY=B得
3.5000
2.6000
3.5556
1.0000
五、算法比较:
Gauss消去法是针对一般的线性方程组,与线性代数中的初等变换解线性方程组方法类似。
高斯消元法的算法复杂度是O(n3);
这就是说,如果系数矩阵的是n×
n,那么高斯消元法所需要的计算量大约与n3成比例。
高斯消元法可用在任何域中。
高斯消元法对于一些矩阵来说是稳定的。
对于普遍的矩阵来说,高斯消元法通常也是稳定的。
LU分解法比较简便迅速,当解多个系数矩阵为A的线性方程组时,LU分解法就显得特别优越,只要对系数矩阵做一次LU分解,以后只要解三角形方程即可。
也可以更具系数矩阵的形状来设计算法。
追赶法只是针对系数矩阵为三对角阵的方程组,因此是一种特殊的方程组。
此方法效率较高,不过不适用于一般的线性方程组。